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ABSTRACT 


MagnetoteUuric  prospecting  is  a  method  of  geophysical  exploration 
that  makes  use  of  the  fluctuations  in  the  natural  electric  and  magnetic  fields 
that  surround  the  earth.  These  fields  can  be  measured  at  the  surface  of  the 
earth  and  they  are  related  to  each  other  by  a  surface  impedance  that  is  a 
function  of  the  conductivity  structure  of  the  earth's  substrata. 

This  report  describes  some  new  methods  for  analyzing  and  interpreting 
magnetotelluric  data.  A  discussion  is  given  of  the  forms  of  the  surface 
impedance  for  various  classes  of  models,  including  one,  two  and  three 
dimensional  models.  Here,  an  n  dimensional  model  is  one  in  which  the 
parameters  describing  the  model  are  functions  of  at  most  n  space  coordinates. 
Methods  are  discussed  for  estimating  the  strike  direction  for  data  that  is 
at  least  approximately  two  dimensional.  A  new  linearized  approach  to  the 
one  dimensional  problem  is  discussed.  Subject  to  the  approximations  of  the 
linearization,  it  is  shown  that  under  the  appropriate  transformations  of  the 
frequency  and  depth  scales,  the  reciprocal  of  the  surface  impedance  as  a 
function  of  frequency  is  equal  to  the  square  root  of  the  conductivity  as  a 
function  of  depth  convolved  with  a  linear  response  function  that  is  somewhat 
like  a  low  pass  filter. 

Included  in  this  report  is  a  comparison  of  several  methods  of  esti¬ 
mating  the  auto  and  cross  power  density  spectra  of  measured  field  data,  and 
of  several  methods  for  estimating  the  surface  impedance  from  these  spectra. 
The  effects  of  noise  upon  these  estimates  are  considered  in  some  detail. 
Special  emphasis  is  given  to  several  types  of  artificial  noise  including 
aliasing,  round  off  or  digitizer  noise,  and  truncation  effects.  Truncation 
effects  are  of  the  most  interest  since  they  depend  upon  the  particular  window 
used  in  the  spectral  analysis. 
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I.  INTRODUCTION 


Magnetotelluric  prospecting  is  a  relatively  new  method  of  “geophysical 
exploration,  although  the  electric  and  magnetic  fields  that  it  employs  have 
long  been  observed.  More  than  a  century  ago  it  was  recognized  by  several 
investigators  that  a  correlation  existed  between  the  variations  in  the  telluric 
currents  and  the  geomagnetic  field.  In  1940  Chapman  and  Bartels  reviewed 
the  various  theories  on  the  relationship  between  these  fields.  In  the  late 
1940's  and  early  1950‘s  several  investigators  such  as  Tikhonov  in  the  USSR; 
Kato,  Kikuchi,  Rikltake,  and  Yokoto  in  Japan;  and  Cagniard  in  France  began 
to  recognize  the  electromagnetic  nature  of  these  fields. 

In  1953  Cagniard  published  a  paper  in  which  he  gave  a  quantitative 
description  of  the  relationship  between  the  electric  and  magnetic  fields  at 
the  surface  of  a  horizontally  layered  earth.  Soon  thereafter  many  people 
began  making  theoretical  and  experimental  contributions  to  the  field  of  mag- 
netotellurics.  By  the  late  1950's,  it  was  recognized  by  several  investigators 
that  the  scalar  impedance  described  by  Cagniard  was  not  sufficient  to  describe 
many  of  the  frequently  encountered  geologic  situations.  For  an  anisotropic 
or  laterally  inhomogeneous  earth,  the  Impedance  becomes  a  tensor  quantity 
(Neves,  1957),  (Rankin,  1960),  (Cantwell,  i960)  ,(Kovtun,  1961) ,  (Rokityanskii, 
1961),  (d'Erceville  and  Kunetz,  1962),  (Bostick  and  Smith,  1962) ,  (Srivastava , 
1963).  Principal  contributors  to  the  growing  body  of  literature  on  magneto- 
tellurics,  in  addition  to  those  previously  mentioned,  include  Berdichevskii, 
Vladimirov,  and  Kolmakov  in  the  USSR,  Porstendorfer  in  Germany,  Adam  and 
Vero  of  the  Hungarian  Academy  of  Sciences,  Fournier  in  France,  and  many 
people  in  the  U.S.A.  and  Canada.  Hugo  Fournier  (1966)  has  a  comprehensive 
history  and  bibliography  of  the  science  of  magnetotellurics. 

The  tensor  relationship  between  the  E  (electric)  and  H  (magnetic)  fields 
at  any  given  frequency  can  be  expressed  as 


1 


(1.1) 


•  — 

— 

— 

r 

X 

Z 

XX 

Z 

xy 

H 

X 

E 

Z 

Z 

H 

_  Yj 

yx 

yy 

y  _ 

where  rectangular  cartesian  coordinates  have  been  Indicated.  This  tensor 
Impedance  Z,  a  function  of  frequency  and  space  coordinates,  depends  upon 
the  conductivity  of  the  earth  In  the  surrounding  area,  and  If  the  horizontal 
wavelengths  of  the  Incident  fields  are  sufficiently  long,  Z  will  be  Independent 
of  time  and  source  polarization.  Therefore,  Z  can  be  a  useful  measure  of  the 
conductivity  structure  of  the  earth,  and  In  fact  It  can  sometimes  be  Interpreted 
almost  completely  In  terms  of  a  simplified  earth  model. 

The  magnetotellurlc  problem  can  conveniently  be  divided  Into  three 
parts:  data  acquisition,  analysis,  and  modeling.  Data  acquisition  Includes 
the  Instrumentation  arid allof  the  field  work  Involved  with  recording  the 
electric  and  magnetic  field  variations.  Analysis  Includes  processing  the 
field  measurements  to  determine  estimates  of  the  Z  tensor  and  other  related 
parameters.  Modeling  consists  of  Interpreting  this  Impedance  tensor  In  terms 
of  a  particular  earth  model. 

The  research  that  went  Into  this  thesis  was  aimed  at  developing  better 
methods  of  magnetotelluric  analysis  and  Interpretation.  The  thesis  itself 
provides  for  the  first  time  a  unified  treatment  of  the  techniques  developed  as 
a  result  of  this  research.  The  treatment  Is  facilitated  by  first  considering 
the  forms  of  theoretical  Impedance  tensors  for  several  classes  of  models. 
Next,  various  methods  are  presented  for  estimating  actual  Impedance  tensors 
from  measured  field  data.  Finally,  the  effects  that  various  types  of  noise 
have  upon  the  impedance  estimates  are  considered. 
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II.  ONE  DIMENSIONAL  MODELS 


In  this  chapter  several  one  dimensional  models  ,  that  is ,  models  which 
have  medium  parameters  that  are  functions  of  only  one  space  coordinate  will 
be  considered. 


A.  Homogeneous  Half  Space  Model 

The  simplest  of  all  possible  models  is  one  in  which  the  earth  is 
considered  to  be  a  homogeneous  ,  isotropic  half  space  of  conductivity  a, 
permittivity  e,  permeability  d  .  Within  any  medium  of  constant  a,  e,  and  U, 
if  we  assume  time  variations  of  the  form  eJU)t,  Maxwell's  equations 


v  x  E  =  -  joja  H 

(2.1) 

v  x  H  =  (a  +  jcue)  E 

— * 

(2.2) 

v  •  H  =  0 

— * 

(2.3) 

v  •  E  =  0 

(2.4) 

combine  to  give  the  vector  Helmholtz  equations 


V 


H 


where 


Y 


2 


2 

=  j'Jiurr  -u )  lie 


(2.5) 


In  rectangular  cartesian  coordinates ,  this  vector  equation  sepa¬ 
rates  ,  so  that  each  of  the  components  of  the  E  and  H  fields  satisfies  the 
scalar  Helmholtz  equation.  Elementary  solutions  to  this  equation  are  of  the 
form 

-(Y  x  +  Y  y  +  Y  z) 
x  v  z 

A  L  * 


where 


3 


(2.6) 


=  jama  -  u)  M-e 


The  general  solution  is  obtained  by  summing  various  elemencary  solutions  with 

different  values  of  A,  Y  ,  Y  ,  and  Y  ,  subject  to  the  constraints  of  equation 

x  y  z 

(2.6) .  Returning  for  the  moment  to  an  elementary  solution,  if  the  coordinate 
axes  are  aligned  such  that  positive  z  is  down,  and  the  direction  of  propaga¬ 
tion  is  in  the  x-z  plane ,  then  the  elementary  solution  is  of  the  form 


*V  ~  ^ zZ  2  2  2 

Ae  2  ;  Y  +Y  =Y^ 

x  z 


(2.7) 


Thus  ,  for  a  homogeneous  plane  electromagnetic  wave  with  its 
direction  oi'  propagation  in  the  x-z  plane ,  each  of  the  comoonents  of  E  and 
H  will  be  of  the  form  shown  in  (2.7) . 

Since  any  homogeneous  plane  wave  can  be  separated  into  TE 
(horizontal  E  field  only)  and  TM  (horizontal  H  field  only)  modes ,  and  since 
the  equations  are  linear  with  respect  to  the  fields,  one  can  consider  the  two 
modes  separately. 

For  the  TE  mode , 


E  =  E  =  0 
x  z 

and  equation  (2.1)  becomes 
_dE  _dE 

-i  ~L+  k  jr*-  =  -  jum(iH  +  jH  +  kH  ) 
oz  ox  x  y  z 

Thus 


Y  E 
z  y 


=  -  jumH 

5C 


-Y 


E 

x  y 


-  jumH 
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H  =  0 
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In  particular, 
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Y  =  jujua 


(2.11) 


Continuity  of  the  tangential  fields  at  the  surface  z  =  0  requires 


V  =  Y 

x  air  x  earth 


(2.12) 


For  a  plane  wave  striking  the  earth  at  a  real  angle  9  measured 
from  normal  incidence , 


Y  .  =  ju)  Vue  sin  9 
x  air 


(2.13) 


Equation  (2.10),  (2.12)  and  (2.13)  together  imply  that 
2  2 

|Y  |  «  l Y  |  .  Therefore  one  may  take 

X 

Y  =  Y  =  /ju)Ua 
z 

Under  these  conditions,  equations  (2.8)  and  (2.9)  give 


(2.14) 


ZTE  ^TM 


(2.15) 


This  implies  that  the  impedance  is  independent  of  the  polarization  of  the 
elementary  solution.  Thus,  any  general  solution  made  up  .of  .elementary  solu¬ 
tions  satisfying  the  conditions  of  equation  (2.14)  will  give  a  scalar  impedance 


(2.16) 


which  will  relate  any  horizontal  component  of  the  total  H  field  to  the  orthogo¬ 
nal  horizontal  component  of  the  E  field . 

Actually  it  is  net  necessary  to  restrict  the  general  solution  for 
the  incident  fields  to  modes  corresponding  to  real  angles  of  incidence,  as 


6 


2  2 

indicated  by  equation  (2.13).  Elementary  solutions  for  which  |y  \  >  ti)  i-*6  will 

X 

still  give  rise  to  total  fields  which  satisfy  equation  (2.16)  provided 

|Y  | 2 «  id ua  (2.17) 

It  is  convenient  to  define  a  parameter  6  ,  called  skin  depth,  for 
conductive  materials  by 

6  =  V  2/ujUa  (2.18) 

Then 

“V  -(1  +  j)z/5  (2.19) 

e  =  e 

Thus  6  is  a  measure  of  the  depth  that  an  electromagnetic  field 
will  penetrate  into  a  conductive  medium .  It  is  the  depth  at  which  the  field 
will  have  been  attenuated  to  1/e  of  its  surface  value . 

If  one  then  defines  the  horizontal  wavelength  X  of  an  elementary 

solution  by 


x 

then  a  statement  that  is  equivalent  to  equation  (2.17)  is  that 
X  »  5 


(2.20) 


In  other  words  ,  the  horizontal  wavelength  is  long  compared  to  the  skin  depth. 

In  summary  then,  for  an  earth  model  consisting  of  a  homogeneous 
half  space  of  conductivity  a,  with  incident  fields  having  horizontal  wave¬ 
lengths  long  compared  to  a  skin  depth,  the  surface  impedance  will  be  given 
by  equation  (2.16). 
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B.  Horizontally  Layered  Model 

The  next  model  that  might  be  considered  is  one  in  which  the  earth 
is  represented  by  a  sc.  of  horizontal  layers,  each  with  a  different  conductivity. 
This  is  usually  known  as  the  Cagniard  model  since  it  is  the  one  that  he  consid¬ 
ered  in  his  classic  paper.  One  assumes  N  layers,  as  shown  in  figure  1,  and 
assumes  elementary  solutions  in  each  layer  of  the  form 


where  d.  is  the  thickness  of  the  ith  layer,  R.  is  a  reflection  coefficient  defined 
1  1 

by 


20i  +  ZR1 


i  -  1,  2,  .  .  .  N-l  (2.23) 


and  ZQi  is  the  characteristic  impedance  of  the  ith  layer.  As  with  the  homo¬ 
geneous  half  space,  the  characteristic  impedance  for  the  TE  mode  is 
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Z1*  =  jtoU/Y 
01 


(2.24) 


zi 

and  for  the  TM  mode  is 

Z™  =  jtud  Y  ,/y2  (2.25) 

01  zi  1 

Again  if  one  assumes  that  the  horizontal  wavelengths  of  the  incident 
fields  are  long  compared  to  a  skin  depth  in  each  layer,  the  two  modes  become 
equivalent  and 


Z  .  =  V joju/a.  (2.26) 

Ol  1 

In  either  case,  one  may  start  at  the  bottom  layer  and  work  up, 
computing  and  Z^  using  the  recursion  equations  (2.23)  and  (2.22)  until  Z^, 
the  surface  impedance,  is  obtained. 

Recall  from  equation  (2. 16)  that  for  a  homogeneous  half  space 


U)H 


Correspondingly,  for  a  layered  model,  it  is  customary  to  define 
an  apparent  conductivity  a  0«)  or  apparent  resistivity  P  (uj)  by 

3  cl 


O  (id) 

a 


urn 

|  Zx(u>)  |2 


_l _ 

p_  (u>) 

a 


(2.27) 


Some  sample  curves  of  apparent  resistivity  versus  frequency  for 
several  models  are  plotted  in  figures  2  and  3.  For  high  frequencies,  0  =° . 

cl  i 

and  for  low  frequencies  ,  Qualitatively  it  appears  that  aa('J))  is  a 

"smoothed  out"  version  of  p(z)  with  frequency  w  being  inversely  related  to 
depth  z. 
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Although  it  is  a  simple  matter  to  obtain  the  surface  impedance 
Zj(uj)  in  terms  of  a{ z)  for  any  layered  model,  the  inverse  problem  of  finding 
a(z)  for  a  specified  Z^(uj)  is  not  so  simple.  It  is  a  nonlinear  problem  that 
in  general  can  be  solved  only  by  using  iterative  techniques.  Computer 
programs  are  available  for  least  squares  fitting  Z(jj)  curves  to  N  layer 
models  (Patrick,  1969). 

Since  a  (u>)  is  a  smooth  curve,  one  might  suspect  that  fine 
details  in  a(z)  cannot  be  determined  from  a  (u>) .  This  in  fact  turns  out  to 

cl 

be  the  case;  only  gross  trends  in  a(z)  can  successfully  be  determined 
from  a  (uj)  . 


C.  General  One  Dimensional  Model 

Consider  the  case  where  o(z)  is  a  continuously  varying  function 

of  z  rather  than  being  restricted  to  a  finite  number  of  homogeneour  layers . 

In  this  case,  the  recursion  equations  (2.22)  and  (2.23)  are  replaced  by  a 

differential  equation  for  Z.  There  are  several  ways  to  obtain  the  differential 

equation.  One  way  is  to  combine  the  recursion  equations  (2.22)  and  (2.23) , 

and  let  Az  replace  d^ ,  and  consider  the  limit  as  Az  approaches  zero.  Another 

simpler  method  pointed  out  by  Swift  (1967)  uses  Maxwell's  equations  directly. 

Consider  the  TE  mode  with  E  =  E  =0.  Equations  (2. 1)  and  (2.2)  give 

x  z 

=  i»UHx 
dE 

Hz 


+ 


3H 
_ x 

3z 


Now 
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,  dE  E  SH 

1  _ y  ,  y  _ x 

H  32  „2  Sz 

x  H 

x 


=  -  ^[J^Hx]  +  -f-[aE  +-^f] 


\r  1 

=  -  ja,U  +  -f-  [aEy  -  — 


jum  +  aZTE  -  ju)U  Zje 


3^Te  ^x  2 

"5? - ium  +  (1  -  -£)  o  z\t 


(2.28) 


Similarly  for  the  TM  mode  one  obtains 


->mU  (1-^)+c,4m 


(2.29) 


Again  when  the  incident  fields  have  horizontal  wavelengths  large- 
compared  to  a  skin  depth,  the  TE  and  l’M  modes  become  equivalent  and 


dZ  „2 

=  -  jwu  +oZ 


(2.30) 


This  differential  equation  is  of  course  nonlinear  in  Z;  however, 
if  one  assumes  a  a(z)  profile  such  that  a(z)  =  a  constant  for  z  >  z^,  then 
Z(z^)  =  •J  and  equation  (2.30)  can  be  numerically  integrated  from 

z  =  z^  to  z  =  0  to  obtain  an  expression  for  the  surface  impedance  in  terms  of 
the  conductivity  profile. 

Thus ,  as  with  the  layered  model,  the  forward  going  problem  of  . 
finding  the  surface  impedance  in  terms  of  a  specified  conductivity  profile  is 
relatively  simple.  Again,  the  inverse  problem  of  finding  the  a(z)  profile  which 
produces  a  specified  surface  impedance  must  be  worked  iteratively. 

D.  Linearized  One  Dimensional  Model 

Consider  the  following  simplification  of  the  one  dimensional  prob¬ 
lem.  Assume  that  the  E  field  as  a  function  of  depth  z  has  the  form 

y(z')  dz' 

E  (z)  =  Ae  0  (2.31) 

X 

where  A  is  independent  of  z  and 


Y  (z)  =  ./  jiuH  o  (z) 


(2.32) 


From  Maxwell's  equations 


dH 


__Z  =  _  0£ 

dz  x 


Integrating  with  respect  to  z  and  noting  that  Hv  must  vanish  as  z  -*  <»  ,  one  has 


eo  dH 


S  df  dz  =  -  S  a  Ex 


dz 


or 
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z 

Y(z')  dz' 

CO 

Hy(a>)  -  Hy(o)  =  -  ^  a(z)  Ae  °  dz 

o 


Thus ,  if  one  defines  the  surface  admittance  Y(uj)  as  being  the  reciprocal  of  the 
surface  impedance,  then 


Y(uj)  = 


H^(z  =  o) 

E  (z  =  o) 
x 


or  /  from  equation  (2.32) 

z 

w  -  'J  jiuU  ^  *J a(z’)  dz' 

Y(io)  =  ^  a(z)  e  °  dz 

o 


(2.33) 


(2.34) 


Now  consider  the  following  transformation.  Let 


and 


e 


y  ujU/2 


Noting  that 


_ 

e  da^  =  V  a  (z)  dz 


(2.35) 


(2.36) 


(2.37) 


and 
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2=0 


<*i  = 


2  =  ®  — *  Ctj  = 


equation  (2.34)  becomes 


—00 


CO 


00 

Y(a2)  =  ^ 


(a,-cu) 

/-T-r  -a  +  j)  e  al  , 

k/o((xJ  e  e  da. 


or 


“a2  2 


-ft  j)  e 


-(a0  -  a.) 


Y(a2)  e  =  J  Vofaj)  e"'1  r 


2  1  -(a, -a.) 

e  xdai  (2.38) 


Noting  from  equation  (2.27)  that 
*/  a  (aj)  =  V  u)U  1  Y  (uj)  I 

d 

equation  (2.38)  gives 

/a  (a)  =  1 J  a  (a)  *  g  (a)  |  (2.39) 

cl 

where 

g(a)=/Te-(1  +  f)e"Va  (2.40) 

Thus ,  under  this  simplified  model ,  which  in  effect  neglects  inter¬ 
nal  reflections  in  the  E  field,  the  apparent  conductivity  can  be  obtained  by 
convolving  the  actual  conductivity  profile  with  a  complex  linear  response 


14 


function  in  a-space.  A  plot  of  the  magnitude  of  g(a)  versus  a  is  shown  in 
figure  4.  The  magnitude  of  g(a)  peaks  up  at  a  =  0  and  decays  as  |<xj  increases. 
Also 

g  (a)  da  =  1 

So  g(a) ,  although  it.  is  complex,  is  somewhat  like  the  response  of  a  low  pass 
filter  with  unity  DC  gain.  This  is  consistent  with  the  earlier  observation 
that  a  (to)  is  a  "smoothed  out"  version  of  a(z)  with  tu  inversely  related  to  z. 

cl 

In  practice,  this  simplified  approach  is  probably  not  very  useful 
by  itself  since  the  assumed  form  of  the  E  field  in  equation  (2.31)  is  not  too 
realistic.  Strictly  speaking  it  is  valid  only  if 

^|  «  V  (2)0(2)  (2.41) 

for  all  z.  On  the  other  hand  this  approach  could  be  quite  useful  for  obtaining 
a  first  guess  to  be  used  in  an  iterative  inversion  scheme.  In  particular  if 
one  simply  assumes  that 

n( a)  a:  a  (a)  (2.42) 

a 

then  frequency  and  depth  may  be  related  through  equations  (2.35)  and  (2.36) 
to  give 


c  (z)  2:  cr  (z  (x>)) 
a  a 

where 
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(2.43) 


(2.44) 


Thus ,  an  approximate  depth  scale  may  be  attached  to  the  frequency  scale  for 

an  apparent  conductivity  curve.  Notice  that  for  a  (to)  =  a  ,  a  constant, 

a 

equation  (2.44)  reduces  to  the  standard  skin  depth,  so  one  may  think  of  z  (as) 

a 

as  sort  of  an  integrated  skin  depth.  In  fact,  for  cases  where  equation  (2.41) 
is  satisfied,  z  (ud)  will  be  the  depth  at  which  the  fields  have  decayed  to  1/e 
of  their  surface  value. 


E.  Generalized  Skin  Depth 

As  suggested  by  the  preceding  paragraph,  it  will  be  useful  to 
generalize  the  idea  of  skin  depth  for  an  inhomogeneous  model.  For  a  homo¬ 
geneous  medium,  the  skin  depth  was  defined  to  be  the  depth  at  which  the 
fields  are  attenuated  to  1/e  of  their  surface  values .  For  an  inhomogeneous 
model,  the  fields  of  course  do  not  have  a  simple  exponential  decay;  however, 
if  one  defines  6  (id)  to  be  the  depth  at  which 


o 


then  6  will  be  a  good  measure  of  the  depth  of  penetration  of  the  fields  and 
as  such  it  may  be  taken  as  the  skin  depth. 

In  the  discussion  of  the  horizontally  layered  model,  the  state¬ 
ment  was  made  that  the  incident  fields  could  be  treated  as  normally  incident 
plane  waves  if  the  actual  horizontal  wavelengths  were  long  compared  to  a 
skin  depth  ir.  each  layer  since  for  that  case 


A  less  restrictive  yet  adequate  requirement  is  that 

6  6 

^  Yz  dz  e:  ^  Y  dz 
o  o 


(2.46) 
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where  6  is  defined  by  equation  (2.45).  Clearly  this  condition  will  exist  pro¬ 
vided 


5 

6 

^  Y  dz 

« 

^  Y  dz 

)  x 

o 

o 

(2.47) 


o  2  2 

where  y"  +  =  Y  =  jwHa  .  But  Y  =  jk  =  j2tt/ X.  Thus 

X  Z  XX 


5Yxdz 


2tt6 

X 


(2.48) 


Also,  from  equation  (2.45),  the  definition  of  6, 
6 

K  Y  dz 


JT 


Thus,  the  condition  (2.47)  will  exist  provided 


2tt6 


«  & 


or 


6  «  x//2rr 

So,  if  the  horizontal  wavelengths  are  long  compared  to  the  skin  depth  de¬ 
fined  by  equation  (2.45),  the  incident  fields  may  in  effect  be  treated  as 
normally  incident  plane  waves  . 

In  conclusion  then ,  the  forward  going  one  dimensional  problem 
is  reasonably  simple.  If  the  incident  fields  are  assumed  to  have  horizontal 
wavelengths  long  compared  to  a  skin  depth,  then  any  horizontal  component 
of  the  H  field  is  related  to  the  orthogonal  horizontal  component  of  the  E 
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field  by  a  scalar  impedance  which  is  related  to  the  conductivity  profile.  The 
inverse  problem  of  estimating  the  conductivity  profile  from  a  measured  surface 
impedance,  while  it  is  nonlinear  has  been  worked  with  some  success  using 
iterative  techniques.  The  simple  linearized  model  discussed  here  should  be 
useful  for  providing  a  first  guess  for  such  iterative  solutions. 
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III.  TWO  AND  THREE  DIMENSIONAL  MODELS 


The  scalar  surface  impedance  discussed  in  the  previous  chapter  is 
not  sufficient  to  describe  the  relationship  between  the  horizontal  E  and  H 
fields  for  a  model  that  has  lateral  variations  in  conductivity.  In  this  chaper 
some  general  relationships  will  be  developed  for  two  dimensional  models 
(models  for  which  c  is  a  function  of  two  space  coordinates,  the  vertical  or 
z  coordinate  and  one  horizontal  coordinate,  say  x)  and  for  three  dimen¬ 
sional  models  (models  for  which  a  is  a  function  of  all  three  space  coordi¬ 
nates).  It  will  be  shown  that  for  these  models  the  impedance  must  be 
expressed  as  a  rank  two  tensor  as  was  indicated  in  equation  (1.1). 

A.  ZTE  and  for  Two  Dimensional  Models 

Consider  again  Maxwell’s  equations  as  stated  in  equations  (2.1) 
through  (2.4).  If  one  assumes  that  the  conductivity  a  is  a  function  of  x  and 
z,  equations  (2.1)  through  (2.3)  are  still  applicable;  however,  equation  (2.4) 
must  be  replaced  by 


r»  .  t~  i?\  —  n 

r  i- i/  KJ 


13. 11 


where  once  again  it  is  assumed  that  displacement  currents  in  the  earth  are 
negligible.  If  one  also  assumes  once  again  that  the  horizontal  wavelengths 
of  the  incident  fields  are  long  compared  to  a  skin  depth,  then  in  the  earth, 
everything  is  essentially  uniform  in  the  y  direction  so  that  equations  (2.1) 
through  (2.3)  together  with  equation  (3.1)  in  component  form  become 


=  -jiqjH 
bz  X 


(3.2) 


~t~£l  =  -jcuii  H 
dx  y 


(3.3) 


=  ju,’H  H 
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(3.4) 


(3.5) 


3H 

__ Z  = 
3z 


3H 
_ x 


3z 


9  H 

_ Y. 

3x 


9  H 
_ x 

9x 


+ 


aE 

x 

3H 
_ z 

3x 

aE 

z 

3  H 
_ z 

3z 


aE 

y 


o 


(3.6) 

(3.7) 

(3.8) 


3(c  Ex)  3(oEz) 
3x  3z 


(3.9) 


Observe  that  the  only  field  components  Involved  In  equations  (3.2) ,  (3.4), 

(3.6),  and  (3.8)  are  E  ,  H  ,  and  H  .  Also,  the  only  components  entering 

y  x  z 

equations  (3.3),  (3.5),  (3.7),  and  (3.9)  are  E  ,E  ,  and  H  .  Thus  it  Is 

x  z  y 

apparent  that  the  two  modes  are  decoupled  and  may  be  considered  separately. 

The  mode  Involving  E  ,  H  ,  and  H  is  usually  called  the  TE  or  E  parallel 

y  x  z 

mode  since  the  E  fLeid  Is  horizontal  and  parallel  to  the  strike.  The  mode 

Involving  E  ,  E  ,  and  H  is  called  the  TM  or  E  perpendicular  mode  since 
x  z  y 

the  magnetic  field  is  horizontal  and  the  electric  field  is  perpendicular  to  the  _ 

strike.  The  strike  Is  the  direction  along  which  there  are  no  variations  in 

the  model  parameters,  In  this  case,  the  y  direction. 

Thus,  for  a  two-dimensional  model  two  impedances  are  required 

to  define  the  relationship  between  the  horizontal  components  of  the  E  and  H 

fields:  Z_=-E  /H  and  Z__  =  E  /H  .  Exact  solutions  for  Z__,(to,x)  and 
TE  y  x  TM  x  y  TE 

ZtM^'3^  ln  terms  g(x'z)  are  not  ^actable  analytically  although  a  few 
approximate  cases  have  been  worked  out.  In  general,  solutions  are  obtain¬ 
able  only  by  using  numerical  methods  such  as  finite  differencing  over  a  two 
dimensional  grid.  Computer  programs  are  available  which  Implement  these 
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techniques  (Patrick,  1969).  It  would  appear  that  the  inverse  problem  of 

finding  o(x,z)  in  terms  of  Z  (uu,x)  and  Z  (cu,x)  could  in  principle  be  solved 

1 L  IM 

using  iterative  techniques  similar  to  those  used  for  the  one  dimensional  Inverse 
problem.  It  is  believed  that  such  a  solution  would  be  unique,  although  no 
proof  is  known.  On  the  other  hand,  the  number  of  calculations  involved  for 
grids  large  enough  to  be  of  interest  is  so  great  that  the  problem  seems  to  be 
out  of  the  range  of  present  day  computers.  Nevertheless,  useful  and  instructive 
information  about  two  dimensional  modeling  can  be  obtained  from  solutions  of 
the  forward  going  problem. 


B.  Z  in  a  General  Coordinate  System  for  Two  Dimensional  Models 
As  was  shown  above,  for  a  two  dimensional  model  the  TE  and  TM 
modes  decouple  when  one  of  the  horizontal  coordinates  is  aligned  with  the 
strike.  It  will  now  be  useful  to  obtain  the  relationship  between  the  tangential 
fields  in  a  coordinate  system  in  which  the  horizontal  axes  are  arbitrarily 
oriented. 

Suppose  that  the  x'  -  y'  coordinate  system  as  shown  in  figure  5 
is  aligned  with  the  strike,  so  that 


and 


ztmh; 


(3.10) 


Suppose  that  the 
respect  to  the  x' 


=  -zteh; 

x-y  coordinate  system  is  oriented  at  an  angle 
-  y'  system  as  shown  in  figure  5.  Then 


9 


(3.11) 


with 


E  = 

E'  cos  8  +  E'  sin  0 

(3.12) 

X 

x  y 

T1  — 

La  — 

-E '  sin  8  +  E'  cos  8 

(3.13) 

y 

x  y 

H  = 

H'  cos  8  +  H‘  sin  8 

(3.14) 

X 

x  y 

H  =- 

y 

-H'  sin  0  +  H'  cos  0 
x  y 

(3.15) 
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or  alternately 


(3.16) 


H'  =  H  cos  6  -  H  sin  9 
xx  y 


H'  =  H  sin  9  +  H  cos  9 
y  x  y 


(3. 


Combining  these  equations  gives 


Thus  If  one  defines 


then 


E  =  E'  cos  9  +  E'  sin  9 
xx  y 

"  (ZTMHpcoS6  +  (-ZTEHy  sine 

=  sin9+H  cos  9)  cos  9  -  Z„„(H  cos  9 

TM  x  y  TE  x 

-  H  s  In  9)  s  in  9 

y 

“  HxC(ZTM  -  ZTE>  sln  9  008  6]  +  Hy  tZTM  cosZ  9  + 
+  ZTEsln2  9] 


E  =  Z  H  +  Z  H 
x  xx  x  xy  y 


and 


Zxx=(ZTM-ZTE)slnec0s9 

=  (^c!S5)sln2e 


Z  =  Z™. .  cos2  9  +  Z  sin2  9 
xy  TM  TE 


f  ZTM +  ZTE  N  .  r  ZTM  ZTE"\  „  n 

=  C - 2 - )  +  C - 2 - >OS2Q 


Similarly  for  the  other  components,  one  obtains 


yZ TM  +  ZTE-n  .  /-ZTM  ZTE  \ 

Zyx  =  -  C  2 )  +  C  2  )  008  2  9 


and 


17) 
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In  summary 


sin  2  0 


In  general,  then,  for  a  two  dimensional  model,  the  tangential 
components  of  E  and  H  are  related  by  a  rank  two  tensor  impedance.  The 
diagonal  terms  of  the  Z  matrix  are  in  general  negatives  of  each  other  and  they 
reduce  to  zero  when  the  axes  are  aligned  with  the  strike. 


C.  General  Form  of  Z  for  Three  Dimensional  Models 

For  three  dimensional  models  where  a  is  a  function  of  all  three 
space  coordinates,  the  six  field  components  are  in  general  all  coupled  to  each 
other,  so  it  is  not  possible  to  separate  the  analysis  into  two  distinct  modes 
as  was  done  for  the  two  dimensional  case.  Nevertheless,  it  is  possible  to 
make  some  general  statements  about  the  relationship  between  the  tangential 
components  of  the  E  and  H  fields. 

It  will  now  be  shown  that  a  rank  two  tensor  impedance  of  the  form 
shown  in  equation  (1.1)  is  unique  and  stable,  subject  once  again  of  course  to 
the  assumption  that  the  horizontal  wavelengths  of  the  incident  fields  are  long 
compared  to  a  skin  depth  in  the  earth.  Also  it  will  be  useful  for  later  purposes 
to  establish  that  in  general  the  vertical  magnetic  field  can  be  expressed  as  a 
linear  combination  of  the  two  horizontal  H  field  components.  That  is, 

H  -•  r  H  +  r  H 
z  zx  x  zy  y 
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(3.21) 


where  r  and  r  are  dimensionless  constants,  subject  also  of  course  to 
zx  zy 

the  assumption  that  the  incident  fields  have  horizontal  wavelengths  long 
compared  to  a  skin  depth.  This  assumption  implies  that  the  incident  fields 
may  be  treated  as  normally  Incident  plane  waves.  This  being  the  case,  the 
incident  fields  can  be  separated  into  two  orthogonal  linearly  polarized  plane 
waves.  Clearly,  for  a  linearly  polarized  normally  incident  plane  wave,  each 
of  the  components  of  the  total  E  and  H  fields  will  be  proportional  to  the  ampli¬ 
tude  of  the  incident  wave.  Thus,  if  the  incident  E  is  linearly  polarized  in  the  x 
direction  then 


E  =  a  E 
x  1  xl 

E  =  a.  E  . 
y  2  xi 

H  =  b  E 
x  1  xi 

H  =  b.E  . 

y  2  xi 


H  =  c  E  . 
z  1  xi 


where  E  is  the  incident  field.  Similarly,  if  the  incident  E  is  linearly  polar¬ 
ized  in  the  y  direction 


a  E 

3  yl 

a  .  E  . 

4  yi 


b_  E  . 

3  yi 

b.E  . 

4  yi 


E  . 
yi 


Since  all  of  the  field  equations  are  linear  with  respect  to  E  and  H,  superposition 
must  hold.  Thus  for  a  general  normally  incident  plane  wave 


E  = 
x 


E  .  +  a„ 
xi  3 


E  . 
yi 
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H  =  [C][B] 

H 

L  yj 

Thus 

[Z]  =  [A]  [B]" 1 
and 

fr]  =  £rzx'rzy]  =  C°3P3' 


(3.22) 

(3.23) 

(3.24) 


So  [Z]  and  [r]  are  defined,  and  E  ,  E  and  H  can  be  expressed  as  Linear 
L  L  x  y  z 

combinations  of  and  H^.  The  onLy  problem  that  might  arise  would  be  if 
[B]  were  singular.  Singularity  of  [B]  implies  that 


b1b4  =  b2b3  (3.25) 

Now  for  any  reasonable  earth  modeL,  the  reflection  coefficient  for  the  mag¬ 
netic  field  at  the  surface  is  almost  unity  so  that  the  total  H  field  is  close  to 
twice  the  incident  field.  Thus,  a  normally  incident  plane  wave  with  E  Linearly 
polarized  In  the  x  direction  (and  hence  H  linearly  polarized  in  the  y  direction) 
will  give  rise  to  total  fields  such  that  will  be  considerably  greater  than  H^. 
Thus 

|b2l  »  |bj| 

Similarly,  for  a  normally  incident  plane  wave  with  E  linearly  polarized  in  the 

y  direction,  H  will  be  somewhat  greater  than  H  .  Thus 

x  y 

|b3l»|b4| 

So  clearly 

I  b2b3 1  >>  lblb4 1 


Comparing  this  with  equation  (3.25)  indicates  that  for  any  reasonable  earth 
model,  [B]  will  not  be  singular,  and  hence  Z  Is  defined  by  equation  (3.23). 

Next  it  will  be  useful  to  observe  the  behavior  of  the  elements  of 
Z  as  the  coordinate  system  Is  rotated.  As  with  the  two  dimensional  model, 
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the  elements  of  Z  in  the  x-y  coordinate  system  will  be  expressed  in  terms 
of  the  elements  of  Z  in  the  x'-y'  coordinate  system  as  shown  in  figure  5. 
The  derivation  is  the  same  as  for  the  two  dimensional  case  except  that 
equations  (3.10)  and  (3.11)  are  replaced  by 

i 


E'  =  Z*  H'  +  Z'  H* 
x  xx  x  xy  y 


and 


(3.26) 


Thus 


E'  =  Z*  H'  +  Z'  H' 
y  yx  X  yy  y 


(3.27) 


E  =  E'  cos  8  +  E'  sin  6 
xx  y 


=  (Z*  H'  +Z'  H* )  cos  6  -I-  (Z1  H’  +Z1  H'  )  sin  8 
xx  x  xy  y  yx  x  yy  y 

=  (Z1  cos8+Z'  sin  8)  H'  +  (Z1  cos  8  +Z1  ‘  sih  0)  H' 
xx  yx  x  xy  yy  y 

=  (Z1  cos8+Z'  sin8)(H  cos8-H  sin  8) 
xx  yx  x  y 

+  (Z'  cos8+Z'  sln6)(H  sin8+H  cos  8) 
xy  yy  x  y 

=  [Z'  cos28+Z'  sln28  +  (Z'  +Z’  )sin8cos0]H 
L  xx  yy  yx  xy  x 

+  [Z’  cos28-Z'  sin28  +  (Z'  -Z'  )sin0cos8]H 

*■  w  vx  w  xx  \ 


So  that 


2  2 

Z  =  Z'  cos  9+Z'  sin  8  +  (Z'  +Z*  )  sin  8  cos  8 

xx  xx  yy  xy  yx 


C 


Z'  + 
xx 


Z' 

-ZZ'N 

y 


2’  -  Z‘ 
y  xx 


^)cos2  8+ 


Z*  + 
xy. 


Z1 


sin  2  8 


Similar  expressions  for  Z  ,  Z  ,  and  Z  are  obtained. 

xy  yx  yy 

Z  =  Z  +Z_  cos  2  9  +  Z„  sin  2  6 
xx  1  2  3 

Z  =  Z .  +  Z_  cos  2  8  -  Z„  s  in  2  9 
xy  4  3  2 


The  results  are 

(3.28) 

(3.29) 
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(3.30) 


where 


■=-z. 

+  Z.  cos  2  0 

-  z0 

sin  2  0 

'yx 

4 

3 

2 

=  z, 

-  Z0  cos  2  0 

-  Z0 

sin  2  0 

'yy 

i 

2 

3 

Z* 

+  Z* 

XX 

yy 

2 

Z' 

-  Z* 

XX 

yy 

2 

Z' 

+  Z‘ 

xy 

yx 

2 

Z' 

-  Z* 

XV 

_ ZS 

If  one  further  defines 


(3.31) 


(3.32) 


(3.33) 


(3.34) 

(3.35) 


Z  (0)  =  Z..  oos  2  0  -  Z_  sin 2  0  (3.36) 

o  5  L 

Then  equations  (3.28)  through  (3.31)  become 


;  =  Z  -  Z  (0  4-45°) 

XX  1  o 

(3.37) 

=  Z.  +  Z  (0) 
xy  4  o 

(3 . 38) 

=  -Z.  +  Z  (0) 
yx  4  o 

(3.39) 

=  Z  +  Z  (e+45°) 
yy  i  o 

(3.40) 

The  function  Z  (0)  traces  an  ellipse  in  the  complex  plane  centered 
o 

on  the  origin  as  0  varies  from,  zero  to  180°..  To  show,  that  this. Is  true  take 

ZQ(0)  =  x  +  jy 

where  x  and  y  are  real.  From  equation  (3.36) 
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x  =  Re  [Z^]  cos  2  0  -  Re  [Z?]  sin  2  6 
=  A  cos  (2  0  -  a) 

y  =  Im  [Z^]  cos  2  0  -  Im  [Z^]  sin  2  0 
=  B  cos  (2  0  -  0) 


letting 

gives 

and 

Thus 

and 

So  that 

or 

or 


2  0  -  o'  -  cp 

0  -  O'  =  cp 
Yo 


x  =  A  cos  cp 

y  -  B  cos  (cp  -  cpQ)  =  C  cos  cp  +  D  sin  cp 
x 

cos  cp  =  - 


r~ 


sincp  =  t  J 1  -  (x/A)‘ 


y  =  c[f  ]  +D[±/l-(x/A)2] 


2  2C  ^  CZ  2  ^2  „  x2  . 

y  ~TXy+^  X  =D 

A  A 

2  2  2  2  2  2  2 
A  y  +  (C  +  D  )  x  -  2AC  xy  -  D  A  =  0 


This  is  the  standard  form  for  an  ellipse  centered  on  the  origin.  Thus  Z  (0) 

o 

traces  an  ellipse  in  the  complex  plane  as  0  varies.  Referring  to  equations 
(3.37)  through  (3.40)  one  observes  that  each  of  the  elements  of  Z  then 
traces  an  ellipse  in  the  complex  plane  as  the  measuring  axes  are  rotated. 


29 


D.  Comparison  of  Z  Matrix  for  Two  and  Three  Dimensional  Models 
As  was  shown  in  the  previous  section,  for  a  three  dimensional 
model,  the  elements  of  Z  trace  ellipses  in  the  complex  plane  as  the  measuring 
axes  are  rotated.  From  equation  (3.18)  one  observes  that  for  two  dimensional 
models,  the  theta  dependent  parts  of  the  elements  of  Z  have  fixed  phases. 

Thus  the  ellipses  degenerate  to  straight  lines  for  the  two  dimensional  case. 

Also  one  observes  from  equation  (3.18)  that  the  diagonal  terms  of  the  Z  matrix 
for  the  two  dimensional  case  have  no  constant  term.  Thus  the  straight  line 
representing  the  locii  of  Z^  and  Z  in  the  complex  plane  passes  through  the 
origin.  Figure  6  illustrates  the  general  form  of  the  locii  of  the  elements  of  Z 
in  the  complex  plane  for  the  two  and  three  dimensional  cases, 

At  the  present  time  solutions  for  the  general  three  dimensional 
problem  are  not  available.  For  this  reason,  it  is  usually  desirable  to  find  one 
dimensional  or  two  dimensional  models  that  approximately  fit  measured  data 
which  in  general  is  of  course  three  dimensional.  It  frequently  happens  that, 
over  some  limited  frequency  range,  measured  data  looks  almost  two  dimensional; 
that  is,  the  Z  ellipses  almost  collapse  to  straight  lines  and  the  diagonal  terms 
of  the  Z  matrix  are  almost  negatives  of  each  other.  This  situation  will  occur 
whenever  there  exists  a  horizontal  direction  along  which  the  conductivity  cross 
section  is  nearly  constant  for  a  distance  of  several  skin  depths  Whenever 
this  situation  exists,  it  is  desirable  to  determine  the  approximate  strike  direc¬ 
tion  and  to  estimate  the  corresponding  Z  „  and  Z  for  comparison  with 

1  £/  J.  IVl 

theoretical  Z's  from  two  dimensional  models. 

Several  methods  have  been  proposed  for  estimating  the  principal 
impedance  axes  [Swift,  1967]  all  of  which  converge  to  the  correct  result  when 
the  data  is  actually  two  dimensional.  From  the  point  of  view  of  the  impedance 
ellipses,  the  most  reasonable  way  seems  to  be  to  take 

ZTE  '  ZTM  ■  Z4  ±  Zi  (3-41> 

where  Z^  is  the  semi-major  axis  of  the  ellipse  Zq(6)  as  defined  in  equation  (3 . 36) 
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This  method  yields  the  principal  directions  as  the  values  of  0  which  maximize 
| Zq( 0)  J  .  A  little  algebra  will  show  that  these  values  of  0  are  given  by  the 
equation 


tan  4  0 


2(x2V  y2y3> 

(x22  +  y22)-(x32+y32) 


(3.42) 


where  and  y.  are  the  real  and  imaginary  parts  of  Z.  respectively,  with 
Z^  being  defined  by  equations  (3.33)  and  (3.34).  Incidentally,  this  method 

gives  the  same  result  as  Swift's  method  of  finding  the  angle  0  which  maxi- 

2  2  2  2 
mlzes  (|Z  I  +  [z  I  }  or  minimizes  Mz  I  +  jZ  I  }. 

1  xy1  1  yx1  1  xx1  yy 

Having  thus  obtained  estimates  of  the  principal  axes  of  the 
impedance  matrix  and  the  corresponding  principal  impedance  values,  it  is 
desirable  to  have  some  measure  of  how  two  dimensional  the  data  actually  is. 
To  accomplish  this,  there  are  two  parameters  that  should  be  considered. 
First,  there  is  the  ratio  of  the  constant  terms  in  the  diagonal  and  off  diagonal 
elements  of  the  Z  matrix.  In  other  words,  the  ratio  where  Z^  and  Z 4 

are  defined  in  equations  (3.32)  and  (3.35).  Second,  there  is  the  ratio  of  the 
minor  axis  to  the  major  axis  of  the  Zq(0)  ellipse.  The  magnitudes  of  both  of 
these  ratios  should  be  small  compared  to  unity  in  order  for  the  data  to  fit  a 
two  dimensional  model. 


E.  Use  of  Hz  for  Determining  the  Strike  Direction 

In  the  previous  section,  an  indication  was  given  as  to  how  one 
might  estimate  the  principal  axes  of  a  measured  impedance  matrix  which  is 
approximately  two  dimensional.  However,  no  method  was  given  for  determin¬ 
ing  which  axis  represents  the  strike  direction.  This  matter  can  be  easily 
resolved  in  terms  of  the  vertical  magnetic  field.  Recall  from  equations 

(3.2)  through  (3.9)  that  H  appears  only  in  the  equations  for  the  TE  mode. 

z 

Thus,  with  the  x'-y’  axes  aligned  with  the  strike,  as  in  section  B  of  this 
chapter 
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Ey  ■  -ZTEHx 


■  riEHx 


In  the  x-y  coordinate  system,  at  an  angle  0  from  the  x'-y'  system 


H  =  H1  =  r'  H*  =  r'  (H  cos  0  -  H  sln0) 
z  z  TE  x  TE  x  y 


Hz  =  Hxr;rec°s6"  Vra51"6 


(3 .43) 


Recall  from  equation.  (3.21)  that  for  the  general  three  dimensional  model 


so  that 


H'  =  r'  H'  +  r1  H' 
z  zx  x  zy  y 


H  =  H'  =  r*  H*  +  r1  H' 
z  z  zx  x  zy  y 


=  r'  (H  cosO-H  sin0)+r'  (H  sin0+H  cos0) 
zx  x  y  zy  x  y 


Thus 


=  H  (r‘  cos  0  +  r'  sin0)+H  (r'  cos0-r'  sin0) 
x  zx  zy  y  zy  zx 


r  =  r’  cos  0  +  r'  sin  9 
zx  zx  zy 


(3.44) 


r  =  r'  cos  0  -  r'  sin  0 
zy  zy  zx 


(3.45) 


Comparison  of  equations  (3.44)  and  (3.45)  with  equation  (3.36)  indicates  that 

r  and  r  like  Z  (0)  trace  ellipses  in  the  complex  plane  as  0  varies, 
zx  zy  o 

However,  the  Important  difference  is  that  the  magnitudes  of  r  and  r  have 

zx  zy 

only  one  peak  every  180°  instead  of  every  90°  like  Zq(0)  .  Furthermore,  one 

observes  from  equation  (3.43)  that  In  the  two  dimensional  limit,  the  angle  0 

that  maximizes  r  is  the  strike  direction, 
zx 

Thus,  when  measured  data  is  approximately  two  dimensional,  the 
angle  that  maximizes  |r  |  should  correspond  to  one  of  the  principal  axes  of 

ZX 
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the  Z  matrix,  and  so  it  is  possible  to  estimate  the  approximate  strike  direction 

and  the  corresponding  Z_  and  Z_. .. 

TE  TM 


IV.  METHODS  FOR  ESTIMATING  THE  Z  MATRIX  FROM  MEASURED  DATA 

Now  that  a  considerable  amount  of  attention  has  been  given  to  the 
forms  of  the  Z  matrix  for  various  classes  of  models  and  to  possible  interpre¬ 
tations  of  Z,  it  is  time  to  consider  some  methods  for  estimating  Z  from 
measured  E  and  H  field  data. 


A.  The  General  Problem 
Consider  the  equation 


E  =  Z  H  +  Z  H 
x  xx  x  xy  y 


where  E  ,  H  ,  and  H  may  be  considered  to  be  Fourier  transforms  of  measured 
xx  y 

electric  and  magnetic  field  data.  If  one  has  two  independent  measurements  of 

"noted  by  E^,  Hyl,  Ex2,  Hx2 


E  ,  H  , 

X  X 

and  H  at  a  given  frequency 

and  H  „ 
y2 

respectively, 

then 

Exi 

Hyl 

Ex2 

Hy2 

XX 

Hxl 

H  . 
yl 

Hx2 

Hy2 

and 

Hxl 

Exl 

Hx2 

Ex2 

2 

xy 

Hxl 

Hyi 

rtx2 

“y2 

provided 


H  ,H  .  -  H  _H  £  0 
xl  y2  x2  yl 


(4.1) 
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Equation  (4.1)  simply  states  the  fact  that  the  two  field  measurements  must 
have  different  source  polarizations.  If  the  two  have  the  same  polarization, 
they  are  not  independent. 

Since  any  physical  measurement  of  E  or  H  will  include  some 
noise,  it  is  usually  desirable  to  make  more  than  two  independent  measure¬ 
ments,  and  then  to  use  some  type  of  averaging  that  will  reduce  the  effects 

of  the  noise.  Suppose  one  has  n  measurements  of  E  ,  H  ,  and  H  at  a 

xx  y 

given  frequency.  One  can  then  estimate  Z^  and  Z^  in  the  mean  square 
sense.  That  is  ,  define 


yk  k  k  k  k 

(E  .  -Z  H  .  -Z  H  )(E  -Z  H  .  -Z  H  .) 
xi  xx  xi  xy  yi  xi  xx  xi  xy  yi 

i  =  1 


where  E^.  is  the-  complex  conjugate  of  E^. ,  etc. ,  and  then  find  the  values 

of  Z  and  Z  that  minimize  'if .  Setting  the  derivatives  of  %  with  respect 
xx  xy 

to  the  real  and  imaginary  parts  of  Z  to  zero  yields 

XX 


n  n 

)  E  H  .  =  Z  /  H  .H*.  +  Z 

XI  XI  XX  t-J  XI  XI 


i  =  l 


i  =  1 


n 

-1  * 

/  H  .H  . 
xy  ->  yi  xi 

i  =  1 


(4.2) 


Similarly,  setting  the  derivatives  of  $  with  respect  to  the  real  and  imagi¬ 


nary  parts  of  Z^  to  zero  yields 


n 


n 


n 

r-1 


/  E  ,H  .  =  Z  /  H  ,H  .  +  Z  )  H  .H  . 
^  xi  yi  xx  ^  xi  yx  xy  Lj  yi  yi 


(4.3) 


i  =1 


i  =  1 


i  =  l 


Notice  that  the  summations  represent  auto  and  cross  power  density  spectra. 
Equations  (4.2)  and  (4.3)  may  then  be  solved  s.  multaneously  for  Z  and 

XX 

Z  .  This  solution  will  minimize  the  error  caused  by  noise  on  E  .  It  is 
xy  x 

possible  to  define  other  mean  square  estimates  that  minimize  other  types 
of  noise.  For  example,  if  one  takes 
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1=1  XX 


♦  - 1  ( 


*  * 

Z  E  .  .  Z 

tPh  .)(  -  H*.  --f- 

zxx  yi  2*  M  Z* 

XX  XX 


* 

H  .) 

yi 


the  resulting  solution  will  minimize  the  error  introduced  by  noise  on  H  . 

There  are  four  distinct  equations  that  arise  from  the  various 
mean  square  estimates.  In  terms  of  the  auto  and  cross  power  density 
spectra ,  they  are 


and 


E  E* 
x  x 


Z  H  E*  +  Z  HE* 
xx  x  x  xy  y  x 


(4.4) 


E  E* 
x  y 


Z  H  E*  +  Z  HE* 
xx  x  y  xy  y  y 


(4.5) 


E  H* 
x  x 


Z  H  H*  +  Z  H  H* 
xx  x  x  xy  y  x 


(4.6) 


E  H* 
x  y 


Z  H  H*  +  Z  H  H* 
xx  x  y  xy  y  y 


(4.7) 


Strictly  speaking,  equations  (4.4)  through  (4.7)  are  valid  only 


if  E^E*  ,  E^E*  ,  etc.  represent  the  power  density  spectra  at  a  discrete 
frequency  u;  ,.  In  practice  however,  Z..  are  slowly  varying  functions  of 
frequency,  and  as  such,  E  E* ,  etc.  may  be  taken  as  averages  over  some 

X  X 

finite  bandwidth.  This  is  fortunate  since  it  facilitates  the  estimation  of 
the  power  density  spectra. 


B.  Estimation  of  Power  Density  Spectra 

There  are  a  variety  of  standard  techniques  available  for  esti¬ 
mating  E  E*  ,  E  E*  ,  etc.  ,  the  auto  and  cross  power  density  spectra, 
x.  x  x  y 

several  of  which  will  be  considered  here  I  In  all  the  cases,  it  will  be 
assumed  that  the  field  components  are  given  as  sampled  time  sequences. 
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1.  One  method  that  was  frequently  used  in  the  past  was  that  of 
using  the  auto  and  cross  correlation  functions  of  the  field  components.  This 
method  makes  use  of  the  fact  that  the  Fourier  transform  of  the  auto  correla¬ 
tion  function  of  a  given  signal  is  equal  to  the  power  density  spectrum  of  that 
signal.  Also  the  Fourier  transform  of  the  cross  correlation  function  between 
two  signals  is  equal  to  the  cross  power  density  spectrum  of  the  two  signals. 
Blackman  and  Tukey  (1958)  have  considered  in  detail  the  various  aspects  of 
estimating  correlation  functions  and  the  corresponding  power  density  spectra 
for  sampled  time  sequences.  They  have  given  careful  attention  to  the  spectral 
windows  that  result  from  truncating  the  time  sequences  and  the  correlation 
functions.  Hopkins  (1966)  and  others  have  used  this  method  for  obtaining 

estimates  of  E  E*  ,  E  E*  ,  etc.  in  magnetotelluric  work.  This  method,  when 
x  x  x  y 

compared  with  the  ones  that  will  be  considered  next,  has  several  disadvan¬ 
tages  .  First  it  is  more  time  consuming  on  the  computer  when  many  cross 
spectra  are  needed.  Second,  it  gives  statistically  correct  results  only 
when  the  signals  are  stationary.  Finally,  it  is  more  susceptible  to  error 
from  the  side  lobes  of  the  spectral  window  when  the  spectra  are  not  reason¬ 
ably  flat.  Blackman  and  Tukey  suggest  that  this  third  disadvantage  can  be 
circumvented  to  some  extent  by  digitally  prewhitening  the  time  sequences 
prior  to  computing  the  correlation  functions . 

2.  Another  method  for  estimating  the  power  density  spectra  of 
the  field  components  begins  by  subdividing  each  of  the  time  sequences  into 
several  blocks .  For  each  data  block  one  computes  the  Fourier  transforma¬ 
tion  to  obtain  estimates  of  E^(u;) ,  E^(u)),  etc.  Then  one  forms  the  products 

EvE*,  EvE*  ,  etc.  Finally,  for  each  frequency,  one  averages  the  products 
xx  x  y 

over  the  several  time  blocks ,  thus  obtaining  time  averaged  estimates  of 

E  E*  ,  E  E*  ,  E  H*  ,  etc.  This  method  is  particularly  well  suited  to  snail 
x  x  x  y  x  x 

digital  computers  since  only  one  time  block  of  data  needs  to  be  stored  in 
memory  at  any  given  time,  and  the  blocks  may  be  quite  small  compared  to  the 
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total  time  sequences.  Also  this  method  is  especially  useful  for  situations 
where  the  signals  contain  noise  bursts  that  are  isolated  in  time .  Such 
noise  bursts  may  arise  from  tape  drop-out,  system  saturation  caused  by 
large  amplitude  signals,  or  many  other  sources.  Such  noise  bursts  are 
often  readily  detectable  so  that  data  blocks  containing  them  may  simply 
be  omitted  from  the  time  average . 

3.  Another  method  for  estimating  the  power  spectra  that  is 
very  similar  to  the  previous  one  consists  of  feeding  the  original  time 
sequences  into  a  bank  of  narrow  band  digital  recursive  filters  spanning  the 
desired  frequency  range.  The  outputs  of  these  filters  are  then  treated  the 
same  as  the  outputs  of  the  block  Fourier  transforms  of  the  previous  method. 
This  recursive  filter  method  has  essentially  the  same  advantages  as  the 
previous  method  together  with  the  additional  advantage  that  it  lends  itself 
quite  readily  to  obtaining  spectral  estimates  equally  spaced  on  a  log  fre¬ 
quency  scale.  This  is  because  the  recursive  filters  may  be  designed  such 
that  they  all  have  the  same  Q  and  have  the  appropriate  spacing  on  the  fre¬ 
quency  scale.  Swift  (1967)  has  used  this  technique. 

4 .  For  the  final  method  to  be  considered  here ,  one  begins  by 

Fourier  transforming  each  of  the  entire  time  sequences  .  The  products 

E  E*  ,  E  E*  ,  etc.  are  then  formed  for  each  harmonic.  Finally  the  products 
x  x  x  y 

are  averaged  over  several  neighboring  harmonics  to  obtain  the  desired  band¬ 
width  .  As  far  as  computation  time  is  concerned ,  this  method  is  quite  effi¬ 
cient  if  one  uses  the  Cooley-Tukey  algorithm  for  fast  Fourier  transforms. 

In  fact,  for  a  given  number  of  multiplications  ,  the  spectral  windows  obtain¬ 
able  by  this  method  are  better  than  those  obtainable  by  any  of  the  Other 
methods  considered  here.  (A  detailed  discussion  of  spectral  windows  is 
included  in  Chapter  V.)  This  method,  like  the  last  one,  lends  itself  readily  to 
constant  Q  estimates  of  the  spectral  density  since  the  number  of  harmonics 
averaged  in  each  band  may  be  made  approximately  proportional  to  the  center 


frequency  of  the  band.  The  primary  disadvantage  of  this  method  is  that, 
compared  to  the  two  previous  methods ,  it  requires  a  fairly  large  number  of 
storage  locations;  in  general  it  requires  a  large  computer.  Actually,  in  a 
modified  form  which  is  not  quite  as  efficient  computationally ,  the  Cooley- 
Tukey  algorithm  is  applicable  to  small  computers.  For  a  detailed  considera¬ 
tion  of  this  algorithm,  see  Cooley  (1965). 

If  then,  by  one  means  or  another,  estimates  of  the  auto  and 
cross  power  density  spectra  are  obtained,  one  can  proceed  to  estimate  the 
elements  of  the  Z  matrix. 


C.  Estimation  of  Z  from  Auto  and  Cross  Power  Density  Spectra 
Consider  equations  (4.4)  through  (4.7).  Under  certain  condi¬ 
tions  ,  these  equations  are  independent  so  that  any  two  of  them  may  be 

solved  simultaneously  for  Z  and  Z  .  Since  there  are  six  possible  distinct 

xx  xy 

pairs  of  equations,  there  are  six  ways  to  estimate  Z  and  Z  .  For  example, 

xx  xy 

the  six  estimates  for  Z  are 


Z 

xy 


Z 

xy 


xy 


(H  E*)(E  E*)  - 
x  x  xy 


(H  E.*)(E  E*) 
x  y  x  x 


(H  E*)(H  E*)  - 
xx  y  y 


(H  E*)(H  E*) 
x  y  y  x 


(H  E*)(E  H*)  - 
xxx  x 


(H  H*)(E  E*) 

X  XXX 


(H  E*)(H  H*)  - 
x  x  y  x 


(H  H*)(H  E*) 
x  x  y  x 


(H  E*)(E  H*)  - 
x  x  x  y 


(H  H*)(E  E*) 
x  y  x  x 


(H  E*)(H  H*) 
x  x  y  y 


(H  H*)(H  E*) 
x  y  y  x 


(H  E*)(E  H*) 
x  y  x  x 


(H  H*)(E  E*) 
x  x  x  y 


(H  E*)(H  H*) 
x  y  y  x 


(H  H*)(H  E*) 
x  x  y  y 


(4.8) 


(4.9) 


(4.10) 


(4.11) 
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(H  E*)(E  H*) 
x  V  x  y 


(4.12) 


and 


xy 


xy 


(H  E*)(H  H*) 

x  Y  y  y 


(H  H*)(E  H*) 

X  X  X  v 


(H  H*)(E  E*) 
x  y  x  y 


(H  H*)(H  F*) 
x  y  y  y 


(H  H*)(E  H*) 
x  y  x  x 


(H  H*)(H  H*) 
x  x  y  y 


(H  H*)(H  H*) 
x  y  y  x 


where  Z  denotes  a  measured  estimate  of  Z 

xy  xy 


(4.13) 


It  turns  out  that  two  of  ther  e  expressions  tend  to  be  relatively  unstable  for  the 

one  dimensional  case,  particularly  when  the  incident  fields  are  unpolarized. 

For  this  case  E  E*  ,  E  H*  ,  E  H*  ,  and  H  H*  tend  toward  zero,  so  that 
xy  xx  yy  xy 

equations  (4.10)  and  (4.11)  become  indeterminant.  The  other  four  expres¬ 
sions  are  quite  stable  and  correctly  predict  Z ^  ~  E for  the  one  dimen¬ 
sional  case,  provided  the  incident  fields  are  not  highly  polarized. 

This  same  thing  is  true  of  the  other  three  impedance  elements 

Z  ,  Z  ,  and  Z  .  In  each  case  there  are  six  ways  to  estimate  Z. . ,  two 
xx  yx  yy  7  ij 

of  which  are  unstable  for  one  dimensional  models  with  unpolarized  incident 
fields .  Also  in  each  case  the  other  four  estimates  are  quite  stable  for  any 
reasonable  earth  model  provided  the  incident  fields  are  not  highly  polarized. 

As  was  mentioned  earlier,  any  physical  measurement  of  E  or  H 
will  necessarily  contain  sonre  noise.  It  is  desirable  now  to  consider  how 
such  noise  will  affect  the  Z  estimates  defined  above.  Suppose  that 


E  = 

E 

+  E 

(4.14) 

X 

xs 

xn 

E  = 

E 

+  E 

(4.15) 

y 

ys 

yn 

H  = 

H 

+  H 

(4.16) 

X 

xs 

xn 

H  = 

H 

+  H 

(4.17) 

y 

ys 

yn 
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where 


and  E  ,  E  ,  H  and  H  are  noise  terms.  If  the  noise  terms  are  all  zero, 
xn  yn  xn  yn 

then  the  four  stable  estimates  of  each  of  the  elements  of  Z  are  the  same,  and 


On  the  other  hand,  when  the  noise  terms  are  nonzero,  the  four  estimates  are 
in  general  different. 

Equation  (4. 13)  for  Z^  corresponds  to  the  one  that  Swift  (1967) 

used.  He  showed  that  his  estimates  of  Z..  were  biased  down  by  random 

noise  on  the  H  signal,  but  were  not  affected  by  random  noise  on  the  E  signal. 

Similar  arguments  for  the  four  stable  estimates  defined  above  indicate  that  in 

each  case ,  two  of  them  are  biased  down  by  random  noise  on  H  and  are  not 

biased  by  random  noise  on  E  (for  example,  equations  (4.12)  and  (4.13)  for 

Z  )  while  the  other  two  are  biased  up  by  random  noise  on  E  and  are  not 
xy 

biased  by  random  noise  on  H  (fo:  example,  equations  (4.8)  and  (4.9)  for 

Z  ).  The  effects  of  the  noise  are  most  easily  seen  for  the  one  dimensional 
xy  _ 

model.  For  this  model,  if  the  incident  fields  are  depolarized  so  that  E  E*  , 
_ _  _  x  y 

E  H*  ,  E  H*  ,  and  H  H*  tend  to  zero,  then  equations  (4.8)  and  (4.9)  for 
_x  x  y  y  x  y 

Z  reduce  to 
xy 


Z  =  E  E*  /  H  E* 
xy  x  x  y  x 


(4.18) 


Equations  (4.12)  and  (4.13)  reduce  to 


Z 

xy 


=  E  H*  /  H  H* 
x  y  y  y 


(4.19) 
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If  one  assumes  that  and  H  are  given  by  equations  (4 . 14)  and 
(4. 17)  and  the  and  Hxn  are  random  and  independent  of  the  signals  and  of 
each  other,  then  the  expected  values  of  the  power  density  spectra  are 


<  E  E*>  = 
x  x 


<  H  H*>  = 

y  y 


<  E  H*>  = 
x  y 


<E  E*  >  +  <E  E*  > 
xs  xs  xn  xn 


<H  H*  >  +  <H  H*  > 
ys  ys  yn  yn 


<H  E*>  =  <E  H*  > 

y  x  xs  ys 


Thus ,  if  the  spectral  estimates  contain  enough  terms  in  the  average  so  that 

the  cross  terms  may  be  neglected  (i.e.  E  E*  ,  etc.  are  negligible),  then 

xs  xn 

equation  (4.18)  gives 


E  E*  +  E  E* 

—  xs  xs  xn  xn  „  „  E  noise  power  v 

Z  =  -  =  Z  (1  +  — — - ,  - ) 

xy  — — - —  xy  E  signal  power 

rl  L  * 

ys  xs 

and  equation  (4.19)  cives 


(4.20) 


Z 

xy 


E  K* 
xs  ys 


H  H*  + 
ys  ys 


H  H* 
yn  yn 


/(1  + 


H  noise  power  . 
H  signal  power  ' 


(4.21) 


Thus  the  estimate  shown  in  equation  (4.20)  is  biased  to  the  high  side  by 
random  noise  on  E  while  the  one  in  equation  (4.21)  is  biased  to  the  low 
side  by  random  noise  on  H.  For  similar  percentages  of  random  noise  on 
E  and  H,  an  average  of  tha  various '.estimates' .hopefully  will  be  better- 
than  any  one  estimate  by  itself.  Also  the  scatter  between  the  various  esti¬ 
mates  should  be  a  good  measure  of  the  amount  of  random  noise  present. 

In  practice  of  course  things  are  not  quite  this  neat  because  the 
assumption  that  the  cross  terms  in  the  average  power  estimates  are  negigible 

may  not  be  valid .  For  example ,  terms  of  the  form  E  H  *  will  not  be 

xn  yn 
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negligible  if  the  two  noises  are  coherent.  Such  might  be  the  case  for  certain 

types  of  instrumentation  noise  or  local  industrial  noise  or  60  cps  power  line 

noise.  Also  terms  of  the  form  E  E*  will  not  be  negligible  if  the  noise  is 

xs  xn 

coherent  with  the  signal  source.  Even  if  all  of  the  noise  terms  are  random 
and  independent  of  the  signals  and  of  each  other,  the  cross  terms  may  not 
be  negligible  if  the  average  power  estimates  do  not  have  enough  degrees  of 
freedom. 
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V.  NOISE  PROBLEMS 


As  was  mentioned  in  the  previous  chapter,  any  physical  measure¬ 
ment  of  E  or  H  will  include  some  noise.  This  noise  may  be  in  the  form  of 
a  constant  bias  caused  by  inaccurate  calibration  of  the  measuring  system, 
or  it  may  be  a  nonlinear  effect  such  as  would  result  from  drift  in  the 
sensitivity  of  the  measuring  system.  On  the  other  hand  many  types  of  noise 
are  Independent  of  the  signal.  These  include  such  things  as  amplifier  noise, 

60  cps  power  line  noise,  digitizer  round  off  noise,  and,  if  the  signals  are 
recorded  in  analog  form,  tape  recorder  noise.  Also,,  there  is  always  the 
possibility  of  having  source  generated  noise.  For  example,  if  the  incident 
fields  include  s^.ne  plane  waves  with  horizontal  wavelengths  short  compared 
to  a  skin  depth  in  the  earth,  the  resulting  surface  fields  may  be  represented 
as  containing  noise. 

In  any  event  one  can  always  represent  the  measured  field  components 
as  sums  of  signals  and  noises  as  indicated  in  equations  (4.14)  through  (4.17). 
The  degree  to  which  the  noise  terms  are  Independent  of  the  signal  terms 
depends  entirely  upon  the  source  of  the  noise.  In  the  cases  where  the  noise 
terms  are  dependent  upon  each  other  or  upon  the  signal  terms,  the  effects  of 
the  noise  upon  the  Z  estimates  vary  according  to  which  estimates  are  used, 
and  according  to  which  signal  and  noise  terms  are  coherent.  No  attempt  has 
been  made  to  catalogue  all  of  the  various  posslb  a  combinations  of  signals  and 
coherent  noises. 

For  the  situation  where  the  noise  terms  are  independent  of  each  other 
and  independent  of  the  signals,  some  interesting  results  can  be  shown. 

A.  General  Incoherent  Noise 

As  was  mentioned  in  the  previous  chapter,  the  various  estimates 
of  the  elements  of  the  Z  matrix  are  biased  either  up  by  random  noise  on  E  or 
down  by  random  noise  on  H.  This  is  caused  by  the  fact  that  the  auto  power 


dent  ity  spectra  are  in  general  biased  up  by  random  noise,  while  the  cross 
powe,-  density  spectra  are  not  biased.  For  example,  suppose  that 


and 


where 


E  =  E  +  E 
x  xs  xn 


H  =  H  +  H 
y  ys  yn 


<  E  E  >  =  0 
xs  xn 


<  H  H  >  =  0 
ys  yn 


<  E  H  >  =  0 
xn  yn 


and  where  the  brackets  <  >  denote  "expected  value  of." 
situation 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 

Clearly,  for  this 


<  E  E*  > 
x  x 


<E  E  >  +  <E  E  > 
xs  xs  xn  xn 


(5.6) 


<  H  H*  > 

y  y 


=  <H  H*  >  +  <  H 
ys  ys 


H*  > 
yn  yn 


(5.7) 


and 


<  E  H  >  =  <  E 
x  y 


H*  > 
xs  ys 


(5.8) 


Equation  (5.8)  Implies  that  the  cross  power  can  be  estimated  to  any  arbitrary 
degree  of  accuracy  by  measuring  the  fields  for  a  long  enough  period  of  time. 
On  the  other  hand,  equations  (5.6)  and  (5.7)  imply  that  the  estimates  of  the 
auto  powers  will  be  biased  regardless  of  the  length  of  time  that  the  fields  are 
measured. 


These  ideas  lead  one  to  consider  an  alternate  approach  to  the 
problem.  Suppose  that  one  performs  two  simultaneous  independent  measure¬ 
ments  of  one  of  the  field  components,  say  E  .  If  the  results  are 

X 
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E  =  E  +  E 
xl  xs  xnl 

(5.9) 

and 

E  =  E  +  E  _ 
x2  xs  xn2 

(5.10) 

where 

<  E  E*  .  >  =  0 
xs  xnl 

(5.11) 

<E  E*  0>  =  0 
xs  xn2 

(5.12) 

and 

< E  ,E*  _>  =  0 
xnl  xn2 

(5.13) 

then 

<E  ,E*  >  =  <E  E*  > 
xl  x2  xs  xs 

(5.14) 

Equation  (5.14)  implies  that  the  E  auto  power  density  spectrum  can  be  estimated 

X 

to  any  arbitrary  degree  of  accuracy  from  two  simultaneous  noisy  measurements 
of  E  if  the  measurements  are  taken  for  a  long  enough  period  of  time  and  if  the 

X 

noises  on  the  two  measurements  are  Independent. 

In  general,  if  one  has  double  measurements  of  either  the  two  tan¬ 
gential  components  of  E  or  the  two  tangential. components  of  H,  one-can  obtain 
estimates  of  the  four  elements  of  the  Z  matrix  that  are  not  biased  by  random  noise, 

B.  Numerical  Noise 

At  this  time  consideration  will  be  given  to  several  specific  types 
of  numerical  noise.  The  term  numerical  noise  as  used  here  refers  to  any  noise 
that  is  artificially  injected  into  the  signal  when  the  latter  is  sampled  for 
numerical  processing. 

1.  Perhaps  the  most  commonly  recognized  form  of  numerical  noise 
is  that  which  is  usually  referred  to  as  aliasing.  In  accordance  with  the  samp¬ 
ling  theorem,  if  a  continuous  function  which  is  sampled  at  a  rate  f  has  any 
frequency  components  greater  than  the  Nyquist  or  folding  frequency  (equal  to 
fQ/2),  these  components  will  be  lost  from  the  sampled  version  of  the  function. 
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If  any  spectral  analysis  Is  performed  on  the  sampled  function,  the  lost  fre¬ 
quency  components  will  appear  folded  down  Into  the  desired  spectrum,  and 
will  of  course  represent  noise.  Since  aliasing  is  a  well  documented  and  well 
understood  phenomenon,  nothing  more  will  be  said  about  It  here  except  to 
note  that  anyone  who  deals  with  magnetotellurlc  data  or  any  other  form  of 
sampled  data  should  be  aware  of  it. 

2.  The  next  type  of  noise  that  will  be  considered  here  is  round  off 
error  on  the  analog- to-digital  converter.  This  type  of  noise  arises  from  the  fact 
that  the  A-D  converter  has  only  a  finite  number  of  discrete  levels.  Typically 
the  signal  passes  through  many  levels  between  sample  points.  For  this  reason, 
the  noise  can  be  characterized  quite  well  as  a  sequence  of  independent  random 
variables  with  amplitudes  ranging  from  e/2  to  -e/2  with  a  flat  distribution  where 
e  is  the  distance  between  adjacent  levels  on  the  A-D  converter.  Thus  the  noise 
spectrum  will  be  flat.  The  total  noise  power  for  2m  data  points  will  be 


e/2 


Total  Noise  Power  =  ^  ^  x2dx  =  eVl2 

-e/2 


(5.15) 


Since  the  spectrum  is  flat,  the  average  noise  power  per  harmonic  would  be 

e2/l  2m  for  m  harmonics.  If  the  signal  spectrum  were  also  flat,  the  average 

2 

signal  power  associated  with  each  harmonic  would  be  about  (Me)  /12m  where 

M  is  the  number  of  digitizer  levels  that  corresponds  to  the  maximum  peak  to 

peak  amplitude  of  the  signal.  Thus,  the  signal  to  noise  ratio  would  be  on  the 
2 

order  of  M  .  In  practice,  it  frequently  happens  that  the  signal  spectrum  Is  not 
flat.  In  this  case  the  expected  signal  to  noise  level  for  a  given  harmonic  is 
about 

^2  Signal  power  In  harmonic _ 

Average  signal  power  per  harmonic 


..2  Siqnal  power  in  harmonic 

_  M  m  ~  ^  ' , — : - ; - 

Total  signal  power 

U  -1  .  .  .  ... 


(5.16) 
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if  there  is  a  total  of  m  harmonics. 

As  an  experimental  check  of  the  digitizer  noise,  the  following 
was  done.  A  typical  set  of  actual  magnetic  field  data  was  selected.  It  was 
Fourier  transformed,  and  each  harmonic  was  multiplied  by  a  theoretical  Z 
computed  from  a  typical  layered  model. .  The  resulting  theoretical  E  was 
Fourier  transformed  back  to  the  time  domain  and  digitized;  that  is,  rounded 
off  to  a  given  number  of  significant  bits.  The  resulting  E  together  with  the 
original  H  were  used  to  compute  an  apparent  resistivity  versus  frequency 
curve.  Figure  7  shows  the  individual  harmonics  of  the  true  E  power  density 
spectrum  along  with  the  expected  digitizer  noise  levels  for  eight  and  twelve 
bit  digitizing.  Figures  8  and  9  show  apparent  resistivity  versus  frequency  for 
the  individual  harmonics  for  eight  and  twelve  bit  digitizing  together  with  the 
true  p  computed  from  the  assumed  model.  Figures  10  through  12  give  the 
corresponding  results  when  the  power  density  spectra  are  first  averaged  in 
bands  of  constant  Q.  From  these  figures,  it  is  seen  that,  as  expected,  the 
apparent  resistivities  computed  from  the  individual  harmonics  have  random 
scatter  when  the  signal  power  is  not  sufficiently  large  compared  to  the  noise. 
Also  as  expected,  the  apparent  resistivities  computed  from  the  averaged  power 
estimates  are  biased  to  the  high  side  by  random  digitizer  noise  on  E.  These 
experimental  results  are  consistent  with  the  theoretical  discussion  of  digitizer 
noise . 

3.  The  next  type  of  numerical  noise  that  will  be  considered  is 
that  which  results  from  truncating  the  time  series  to  a  finite  length  T.  Suppose 
that  one  of  the  field  components  has  an  amplitude  that  is  described  by  f(t)  fo- 
all  time.  The  Fourier  transform  F(&)  is  then  given  by 

00 

F(w)  =  ^  f(t)  e'jWt  dt  (5.17) 

-00 

It  is  then  desired  to  approximate  F(uu)  by  F  (ntu  ),  a  Fourier  series  representation 
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of  f(t)  over  some  time  interval  T.  Thus 
_  T/2  -jncu  t 

F  (nouj  =  J  f(t)  e  '  0  dt;  n  =  0  1 ,  ±  2,  . . .  (5.18) 

-T/2 


where  ou  =  2n/T. 
o 

Notice  that 


00 

F(u>)  =  J  f(t)  d(t)  e"jWt  dt 


(5.  IS) 


where  d (t)  =1  for  1 1 1  <  T/2  and  d(t)  =  0  for  1 1  j  >  T/2.  Thus,  from  the 
convolution  theorem 


00 

F  (nwQ)  =  ^  F(cu)  D(n(Uo  -  to)  duu 

-  00 

where 


00 


D(uj)  =  ^ 

dft)  e  ^  dt 

-00 

T/2 

-$ 

-jwt  .. 
e  dt 

-T/2 

i 

_  2n 

s  in  (tui/uj  ) 
o 

w 

o 

(nw/w  ) 

0 

(5.20) 


(5.21) 


D((«)  is  usually  called  the  spectral  window  since  the  observed  spectrum  F(w)  is 
equal  to  the  true  spectrum  F(w)  convolved  with  D(u>) . 

The  spectral  window  defined  by  equation  (5.21)  is  actually  not  very 
desirable  since  the  side  lobes  go  off  only  as  l/iu.  A  better  window,  usually 
known  as  the  Hanning  window  is  obtained  by  letting 


d(t)  =  { 


.5  + 
0 


5  cos  u)Qt 


;  1 1  |  <  T/2 
;  1 1  j  >  T/2 


(5.22) 
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Then 


D(u>) 


(.  5  +  .5  cos  u)  t) 
o 


-1/2 

2 

u)  sin  (niu/u)  ) 
o  o 

,2  2  " 

U)  (ou  —  0)  ) 


-jwt 

e 


dt 


(5.23) 


The  main  lobe  of  this  spectral  window  is  twice  as  wide  as  the  main  lobe  of 

3 

the  previous  one;  however,  the  side  lobes  go  off  as  l/u)  .  The  two  windows 
are  compared  in  figure  13. 

One  could  define  windows  that  have  even  smaller  side  lobes; 
however,  they  would  necessarily  have  wider  main  lobes,  and  as  will  be  seen 
later,  this  is  not  desirable.  The  Hanning  window  seems  to  be  an  adequate 
compromise  between  main  lobe  width  and  side  lobe  height. 

One  is  then  faced  with  the  fact  that  any  physical  estimate  of  the 
power  density  at  a  particular  frequency  u>  is  necessarily  a  weighted  average 
of  the  true  power  density  over  a  band  of  frequencies,  the  weighting  function 
being  the  spectral  window  D(uj).  If  the  Impedance  function  that  one  is  attempt¬ 
ing  to  estimate  does  not  change  significantly  over  the  bandwidth  defined  by 
D(u>) ,  then  the  estimate  will  not  be  corrupted  by  the  truncation  effects.  In 
practice,  however,  the  impedance  does  change  some  so  that  there  will  be  some 
truncation  noise.  The  problem  is  particularly  severe  if  the  power  density 
spectra  have  resonant  peaks  or  other  steep  slopes.  If  one  is  attempting  to 
estimate  the  power  density  near  the  bottom  of  a  steep  slope,  the  contributions 
from  the  side  lobes  of  the  spectral  window  may  be-  significant  compared  to  the 
contribution  from  the  main  lobe.  This  effectively  broadens  the  bandwidth  over 
which  the  impedance  function  must  not  change. 

These  considerations  lead  one  to  inquire  into  the  spectral  behavior 
of  the  E  and  H  fields  used  in  magnetoteliuric  surveying.  One  would  hope  that 
the  general  shape  of  the  spectra  of  the  incident  E  and  H  fields  might  be  more 


50 


or  less  independent  of  time  and  space  coordinates.  If  this  were  the  case,  then 
the  measured  total  H  field,  which  is  close  to  twice  the  Incident  H  field,  would 
also  be  reasonably  stationary  with  respect  to  time  and  space  coordinates,  and 
hence  it  could  be  prewhitened.  On  the  other  hand,  the  total  E  field  is  a 
strong  function  of  the  local  conductivity  structure  and  hence,  although  it  would 
be  stationary  with  respect  to  time,  the  shape  of  the  spectrum  would  change 
from  one  location  to  the  next  as  the  conductivity  structure  changes.  But  still, 
the  surface  impedance  Z  is  a  well  behaved  function  of  frequency  and  as  such 
one  would  expect  that  if  the  E  signals  were  passed  through  the  same  filters  as 
were  designed  to  prewhiten  the  H  signals,  the  resulting  filtered  E  signals 
would  have  a  reasonably  well  behaved  spectrum.  This,  in  general,  turns  out 
to  be  the  case.  Actually,  as  is  indicated  by  equations  (2.16)  and  (2.27),  Z 
tends  on  the  average  to  be  proportional  to  the  square  root  of  frequc...  .  so 
that  an  optimum  filter  for  E  would  differ  from  the  H  filter  by  a  factor  of  l//uT. 

With  these  ideas  in  mind,  a  study  was  made  of  the  spectra  of 
some  actual  H  field  data  recorded  in  central  Texas.  Figures  14  and  15  give 
composite  plots  of  and  H  power  density  spectra  obtained  from  104  different 
data  samples  recorded  at  five  different  sites  in  central  Texas.  This  data  was 
recorded  by  D.  It.  Word,  and  a  magnetotelluric  interpretation  of  the  data  is 
given  by  him  (Word,  1969).  From  these  figures  it  is  apparent  that  ac  least  for 
the  locations  and  times  Involved  here,  the  general  shape  of  the  H  power  density 
spectra  is  fairly  well  defined.  However,  there  are  some  definite  resonant  peaks 
(for  example,  around  .07  cps  and  around  2.5  cps)  that  appear  in  some  of  the 
spectra  but  are  absent  from  others.  These  results  are  consistent  with  those 
obtained  by  other  investigators  (Hopkins,  1966),  (Bleil,  1964).  It  is  believed 
that  if  the  analog  H  signals  are  prewhitened  according  to  the  general  trends 
shown  in  figures  14  and  15,  the  Hanning  window  can  be  used  without  encounter¬ 
ing  any  side  lobe  difficulties  except  perhaps  immediately  adjacent  to  the 
observed  resonances. 
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In  order  to  get  an  estimate  of  the  effects  of  truncation  upon  the 
individual  harmonics  of  the  Fourier  spectra,  consider  the  following  problem. 
Assume  a  one  dimensional  case  with 


E(uu)  =  Z(u>)  H(u>)  . 


Suppose  that  the  H  signal  is  prewhitened  with  a  filter  that  has  a  response 


Fh(ou),  and  the  the  E  signal  is  passed  through  a  filter  whose  response  is 
Fe(iu)  which  may  or  may  not  be  the  same  as  F^(uu).  Assume  that  the  outputs 
of  these  filters  are  Hq(uu)  and  Eo(uu)  respectively,  so  that 


Z(uu) 

Then  define 

G  (uj) 


E(uu)  _ 

Eq(u»)  FhM 

H(uu) 

Hq(uu)  FeN 

E0W> 

Z(uu)  F£(uu) 

HoW 

fhw 

(5.24) 


(5.25) 


Thus  G(uu)  is  the  ratio  of  the  prewhitened  electric  and  magnetic  field  signals 
and  will  be  equal  to  Z(iu)  if  the  two  prewhitening  filters  are  the  same.  If  the 
prewhitened  signals  are  then  sampled  and  a  Fourier  series  analysis  is  per¬ 
formed  on  each,  the  results  will  be 


and 


H  (nuu  -  id)  D(uu)  duu 
o  o 


Eo(ni«o-  i d)  D(uu)  duu 


H0(nuuQ-  u))  G(nuuQ-  id)  D(uu)  duu 


(5,26) 


(5.27) 


where  D(id)  Is  again  the  spectral  window  used.  If  G  were  constant  over  the 
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(5.35) 


G  (mi)  )  =  G  (nil)  )  -  G‘  (nil)  ) 
o  o  o 


^  H^niw^  -  iu)  D(w)  wdu) 

-00 

00 

\  H  (mi)  -  tw)  D(<a)  diu 
J  o  o 


and 


error  =  - 


G'  (nil)  ) 
o 

G(nuio) 


^  Ho(nu)o)-  at)  D(io)  uidu) 


-00 


m 

C  H  (nu)  -U))D(uu)  dti) 
J  o  o 


(5.36) 


Assume  now  that  the  Integrals  In  equation  (5.36)  can  be  approximated  by  sum¬ 
mations  of  the  form 


where 


and 


Then 


H(niu  )  =  {  H  (mu  -  ou)  D(iu)  dtu 
o  J  o  o 


m 

=  L 

i  =  -m 


H  iD1 
on, i  i 


H  .=  H  (nil)  -  iAiy) 
on,i  o  o 


m 

y  D 

L  i 

i  =  -m 


=  1 


W  (nu)Q)  -  C  Ho(nu)o- id)  D(iu)  ii)du) 


(5.37) 

(5.38) 

(5.39) 


m 


-l 


H  D  — 
on , i  i  m 


l  =  -m 


(5.40) 
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where  0  is  the  maximum  value  of  ui  contributing  to  the  integral.  That  is, 

«j  =  fl  corresponds  to  i  =  m. 

Assume  that  for  any  fixed  n,  H  .  are  a  sequence  of  2m +  1 

on,  l 

independent  complex  random  variables  whose  components  have  normal  distri¬ 
butions  and  zero  means.  Also  assume  that  the  original  H  signal  was  pre- 
whltened  enough  so  that  over  any  given  band,  the  expected  value  of  the  power 
density  is  constant.  That  is. 


,H*  .>  =  a  ;  1  =  0,  ±1,  ±2,...  ±m 

on ,  l  on ,  i  n 


(5.41) 


Since  H  .  are  assumed  to  be  independent  and  have  zero  means, 
on,  i 

<H  >  =  0 
on,i 


(5.42) 


<Hn  .  Hj  ,‘>  =  0  ;  i#ij 
on,i  on,] 


(5.43) 


Since  H  .  are  complex  normal  random  variabiea,  H  (nio  )  and  W(nuu  ) ,  which 
on , l  o  o 

are  linear  combinations  of  H  .,  must  also  be  complex  normal  random  variables. 

_ on , i  _ 

Consider  the  statistics  of  H  (nuu )  and  W(nu> ). 

o  o 

„  m 

<  H  (nu>  )  >  =  <  Y  H  ,  D.  > 
o  L  on , i  1 

i=-m 


in 

=  .  y  <h  ,>0, 

L  on , i  i 


l  =  -m 


=  0 


(5.44) 


<H(n.o)H*(»o)>=<(£  HonlDl)(|;  H*nJDr)> 

i=-m  j=-m 


(Eq.  cont'd  on  next  page) 
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<  W(no>o)  > 


<W(nou  )W*(nu)  ] 
o  o 


<H  (nuu  )  W*  (nu)  ) 
o  o 


m  m 


l 


D.D,  <H  .  H  > 


L  L,  "i  j  **on,i*‘on,j 
i  =  -m  j  “-m 


m 

=  Y  D.2<H  ,H*  > 

L  i  on,t  on.i 

i  =  -m 


m 

2  r  _  2 


=  a  )  D. 
n  L  i 

i.  =  -m 


(5.45) 
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=  <y  H  D  —  > 

L  on,l  i  m 


i  =  -m 


m 

--  1 
m  Li 


<H  ,  >D  i 
on,  l  i. 


l=-m 


(5.46) 
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>  =  <CI  HoMDi»XI  'Jon,)DJ^> 

j=-m 


l  =  -m 


2  m  m 


■  Cl)  l  l  " 

i  =  -m  j  =  -m 


<H  .H*  .> 

on,i  on,j 


2  m 

(I)  l  ‘V 

1  =  -m 


(5.47) 


m  m 

=  <(I  Hon,lDl)(Z  Hcn,JDj  "m") 


i  =  -m 


j=-m 


(Eq.  cont'd.) 
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m  m 


-i  l  l  «¥>, 

i  =  -m  j  =-m 


<  H  .  H  .  > 
on, l  on,j 


-a  ,*  i  ,D* 


i  =  -m 


I  = 


so  that 


l  =  -m 


<  H  (nuu  )  wj  )  >  =  0 
o  o 


(5.48) 


2  2 

If  the  spectral  window  D(uu)  Is  an  even  function  of  tu,  then  =  D  and 


(5.49) 


(5.50) 


Thus ,  H  (nioQ)  and  W  (n(UQ)  are  Independent  complex  normal  random  variables 
with  zero  means  and  with  variances  defined  by  equations  (5.45)  and  (5.47). 

A  standard  exercise  In  random  variable  theory  Indicates  that  If  two 
random  variables  X  and  Y  are  normal  and  Independent  with  zero  means  and 
equal  variances,  then  the  function 

r* - 

v  ^  Vx2  +  y2 

has  a  Rayleigh  distribution  (see  for  example,  Papoulls,  1965).  Since  H(nu)Q) 
and  W(n<juQ)  are  complex  normal  random  variables  with  zero  means,  their 
real  and  Imaginary  parts  satisfy  the  conditions  of  the  above  exercise.  Thus 
|H(nu)Q)|  and  |W(nuuo)|  must  have  Rayleigh  distributions. 


that  If 


Another  standard  exercise  In  random  variable  theory  Indicates 


V  ^  X/Y 


then  the  distribution  of  V  Is  given  by 


oo  yy  O  » 

Fv(v)  =  ^  J  Ux,y)  dx dy  +  ^  J  f^(x,y)  dxdy  (5.51) 


O  -o# 
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where  f^(x,y)  is  the  joint  density  of  X  and  Y  (see  Papoulis,  1965) .  If  X  and 
Y  are  independent  and  have  Rayleigh  distributions  then  equation  (5.51)  becomes 


r  P 

PV(V)  =  5  $ 


fx(x)  f^(y)  dxdy 


o  o 


=  ^  fy(y)  Fx(yv)  dy 


■s 


JL 

2 
a 

o  y 


-y2/2a2  -(yv)2/2a  2 

e  y  [1  -e  x 


]  dy 


=  1  - 


1 


va  2 
x 


(5.52) 


2  2  2  2 
where  a  =  <x>  and  a  =  <Y  >. 
x  y 

Now  recall  that  |H(nu)Q)|  and  |W(muo)  |  have  Rayleigh  distributions, 
From  equations  (5.45)  and  (5.47) 

,2 


-  m 

2  v1  ^2 


<|H(n»o)|>  =  on  l  D‘ 
i  =  -m 


and 


m 


<|W<»0,|2>=(|)  o2  l  l2  D2 

i  =  -m 

Also,  from  equations  (5.36),  (5.37),  and  (5.40), 


(5.53) 


(5.54) 


error  = 


G(nu)Q)  W(m«o) 

G'lnio  )  H(nw  ) 
o  o 


(5.55) 


Recall  that  F^(v)  is  by  definition  the  probability  that  V  £  v.  Thus,  equations 
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(5.52)  through  (5.55)  combine  to  give 


where 


P  { |  error  |  £  e }  =  1 


1 


1+  (e 


G(nu>o) 

G*(n<jJo) 


m 


i 

l  =  -m 


<*>2  s 

l  =  -m 


12D 


2 

l 


(5.56) 


(5.57) 


and  where  P{X}  denotes  the  probability,  then  X  occurs.  Recall  from  equations 
(5.37)  through  (5.40)  that  the  summations  arose  as  approximations  to  integrals. 
If  one  now  lets  m  -*  •'  and  passes  back  to  the  integral  formulation,  equation. 
(5.57)  becomes  m 

^  D2(uu)  du) 

K2  =  rjr—9  - - -  (5.53) 

^  uu  DZ(uu)  duu 


Since  D(u))  has  been  assumed  to  be  on  even  function  of  uj  ,  and  since 
P  { |error|  >  e  }  =  1  -  P  (|error|  s  e  ) 
one  has  from  equation  (5.56) 


P  { |  error  |  >  e  } 


1  +e 


G(nu>o)  |‘ 
!  G'(muo) ) 


D(u>)  dio 


w 

^  ui2D(uj): 


duu 


(5.59) 


If  D(tu)  is  a  block  window  of  width  BW,  equation  (5.59)  becomes 
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or,  for  the  worst  case 


I  G(nu>  )  j  4nuu 
o  I  o 

;  |  '  T" 


(5.67) 


Thus,  for  this  case  equation  (5.62)  becomes 

1 


P  { |error|  >  e  3  = 


,  .  16  2  2 
1  +  T  n 


(5.68) 


Equation  (5.68)  indicates  that,  as  expected,  the  probability  of  the 
error  being  greater  than  e  goes  down  as  e  increases.  Also,  as  expected,  the 
probable  error  decreases  as  n,  the  harmonic  number.  Increases.  The  latter 
result  is  expected  since  |  AG(uj)/G(uj)  j  should  be  proportional  to  the  percentage 
bandwidth  of  the  window,  which  in  turn  is  inversely  proportional  to  the  harmonic 
number.  The  results  of  equation  (5.68)  are  summarized  in  figure  16. 

It  is  doubtful  that  the  results  shown  in  figure  16  are  useful  quanti¬ 
tatively  because  of  the  assumptions  made  about  the  form  of  Hq(u))  ,  the  pre¬ 
whitened  magnetic  field  signal.  In  particular,  it  was  assumed  that 
2 

ctH  =  <  Hq(ou)  Hq(uj)*>  is  independent  of  frequency.  In  practice  this  is  not 
attainable  since,  as  noted  earlier,  magnetotelluric  signals  are  not  really  sta¬ 
tionary. 

In  order  to  get  some  type  of  estimate  of  the  effects  of  truncation 
upon  realistic  data,  the  following  experiment  was  performed.  A  typical  set  of 
actual  magnetic  field  data  was  selected.  It  was  Fourier  transformed  using  the 
Hanning  window,  and  each  harmonic  was  multiplied  by  a  theoretical  Z  computed 
from  a  typical  layered  model.  The  resulting  theoretical  E  was  Fourier  trans¬ 
formed  back  to  the  time  domain.  This  E,  together  with  the  original  H,  were 
truncated  to  some  fraction  of  the  original  length.  The  resulting  truncated  E 
and  H  signals  were  Fourier  transformed,  and  apparent  resistivities  were  com¬ 
puted  from  each  harmonic.  These  apparent  resistivities  were  then  compared 
with  the  true  apparent  resistivities  for  the  assumed  model.  This  experiment 
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v'as  repeated  for  several  different  models  with  several  different  original  data 
lengths  and  several  different  truncated  data  lengths.  The  results  showed  a 
very  definite  trend.  In  each  case,  the  apparent  resistivities  computed  from 
the  first  eight  to  ten  harmonics  were  significantly  in  error.  The  higher 
harmonics  showed  very  little  error.  The  amount  of  error  on  the  first  few 
harmonics  depended  significantly  upcn  how  white  the  spectra  were;  the 
whiter  spectra  had  less  error.  Figures  17  and  18  show  the  results  of  a 
typical  run.  Figure  17  shows  the  individual  harmonics  of  the  spectrum  of  a 
truncated  H  signal.  Figure  18  shows  the  corresponding  apparent  resistivities 
along  with  the  true  apparent  resistivity  curve  for  the  assumed  model. 

These  considerations  lead  one  to  believe  that  the  first  few  har¬ 
monics  of  a  Fourier  spectrum  of  typical  magneto  telluric  data  are  likely  to  be 
corrupted  considerably  by  truncation  error. 
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VI.  CONCLUSION 


The  methods  of  analysis  discussed  in  the  foregoing  chapters  have  been 
implemented  in  a  digital  computer  program  constructed  by  the  author  for  use  on 
the  CDC  6600  computer  at  The  University  of  Texas  Computation  Center.  This 
program  estimates  the  power  density  spectra  of  sampled  E  and  H  signals  by 
computing  the  Fourier  transforms  of  the  sampled  data  using  the  Cooley-Tukey 
algorithm,  and  averaging  the  resulting  auto  and  cross  powers  in  frequency  bands 
of  constant  percentage  bandwidth  as  discussed  in  Chapter  IV,  Section  B.4.  The 
Hanning  window  discussed  in  Chapter  V,  Section  B.3  is  used  for  the  Fourier 
transforms.  The  elements  of  the  Z  matrix  are  then  estimated  from  the  power 
density  spectra  using  the  techniques  described  in  Chapter  IV,  Section  C.  The 
principal  axes  are  then  determined  in  accordance  with  the  discussion  in 
Chapter  III,  Sections  C  and  D.  Also  the  approximate  strike  direction  is  deter¬ 
mined  from  the  vertical  magnetic  field  as  discussed  in  Chapter  III,  Section  E. 

As  diagnostics,  the  tensor  coherency  mentioned  in  Chapter  IV,  Section  C,  and 
the  two-dimensionality  parameters  mentioned  in  Chapter  III,  Section  D,  are 
determined. 

This  program  has  been  used  extensively  for  analysis  of  magnetotelluric 
data  recorded  in  central  Texas  by  Darrell  Word.  Samples  of  the  results  are 
given  by  Word  (1969).  For  most  of  the  data  analyzed  using  this  program,  the 
resulting  surface  impedance  estimates  have  been  consistent  and  repeatable. 

In  areas  where  the  geology  is  reasonably  one  dimensional,  these  surface  imped¬ 
ances  have  been  successfully  interpreted  in  terms  of  horizontally  layered  models. 
The  resulting  resistivity  profiles  have  agreed  quite  well  with  independent  obser¬ 
vations  such  as  resistivity  well  logs  in  cases  where  the  latter  have  been  avail¬ 
able. 

In  areas  where  the  geology  is  more  complex,  particularly  where  it  is 
highly  three  dimensional,  interpretation  has  not  been  so  successful.  However, 
even  here  the  surface  impedance  estimates  have  been  fairly  repeatable.  Thus, 
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it  is  believed  that  now,  perhaps  for  the  first  time,  the  surface  impedances 
have  been  measured  more  accurately  than  they  can  at  present  be  interpreted. 
For  this  reason,  it  is  believed  that  future  contributions  to  the  science  of 
magnetotellurics  must  come  in  the  area  of  interpretation. 
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Fig.  1  Description  of  N  Layer  Model 
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.Fig.  2  Sample  Two  Layer  Apparent  Resistivity  Curves 
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Frequency  ( cps) 

Fig.  3  Sample  Five  Layer  Apparent  Resistivity  Curves 


Linearized  One  Dimensional  Prob 


Fig.  5  Relative  Orientation  of  x'-y'  and  x-y  Coordinate  Systems 
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Fig.  10  Average  E  Power  Density  Spectrum  for  Digitizer  Noise 
Test  with  Expected  Noise  Level  for  Eight  Bit  Digitizing 


•Mt 


75 


Apparent  Resistivity  ( ohm  -meters  ) 


+  +  from 

overaged 

spectra 

-  from 

assumed 

model 

10  10  10  I' 

Frequency  (  cps  ) 

Fig.  12  Apparent  Resistivity  versus  Frequency  for  Average  Power 
Density  Spectra  with  Twelve  Bit  Digitizing 


Fig.  13  Comparison  of  the  Block  and  Hanning  Spectral  Windows 
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.Magnetotellurlc  prospecting  Is  a  method  of  geophysical  exploration  that  makes  use  of 
the  fluctuations  In  the  natural  electric  and  magnetic  fields  that  surround  the  earth. 
These  fields  can  be  measured  at  the  surface  of  the  earth  and  they  are  related  to  each 
other  by  a  surface  Impedance  that  Is  a  function  of  the  conductivity  structure  of  the 
earth's  substrata. 

This  report  describes  some  new  methods  for  analyzing  and  Interpreting  magnetotellurlc 
data.  A  discussion  Is  given  of  the  forms  of  the  surface  Impedance  for  various  classes 
of  models,  including  one,  two  and  three  dimensional  models.  Here,  an  n  dimensional 
model  Is  one  In  which  the  parameters  describing  the  model  are  functions  of  at  most  n 
space  coordinates.  Methods  are  discussed  for  estimating  the  strike  direction  for  data 
that  Is  at  least  approximately  two  dimensional,  new  linearized  approach  to  the  one 
dimensional  problem  Is  discussed.  Subject  to  the  approximations  of  the  linearization, 
It  Is  shown  that  under  appropriate  transformations  of  the  frequency  and  depth  scales, 
the  reciprocal  of  the  surface  Impedance  as  a  function  of  frequency  Is  equal  to  the 
square  root  of  the  conductivity  as  a  function  of  depth  convolved  with  a  linear  response 
function  that  Is  somewhat  like  a  low  pass  filter. 

Included  In  this  report  is  a  comparison  of  several  methods  of  estimating  the  auto  and 
cross  power  density  spectra  of  measured  field  data,  and  of  several  methods  for 
estimating  the  surface  impedance  from  these  spectra.  The  effects  of  noise  upon  these 
estimates  are  considered  In  some  detail.  Special  emphasis  Is  given  to  several  types 
of  artificial  noise  Including  aliasing,  round  off  or  digitizer  noise,  and  truncation 
effects.  Truncation  effects  are  of  the  most  Interest  since  they  depend  upon  the 
particular  window  used  In  the  spectral  analysis. 
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